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I.  INTRODUCTION 


The  Hydrographic  Sciences  Group  at  the  Nava!  Postgraduate  School  (NPS), 
conducted  its  first  Sea  Floor  Benchmark  Experiment  in  May  1985.  The  experiment 
basically  was  to  deploy  a  benchmark  on  the  ocean  bottom  and  to  establish  its  position 
with  geodetic  accuracy.  This  was  to  be  accomplished  using  GPS  and  acoustic  range 
data  (Kumar,  et  al.,  1984). 

Beyond  the  goal  of  the  Sea  Floor  Benchmark  Experiment  in  May  1985  was  a 
higher  goal  of  establishing  and  continuing  a  GPS  research  base  at  NPS.  In  order  to 
reach  this  goal  it  has  become  evident  that  a  library  of  GPS  programs  should  be 
established  at  NPS  capable  of  reducing  GPS  data,  independent  of  outside  sources,  into 
a  format  compatible  for  analysis. 

The  Texas  Instruments  receiver  (TI-4100)  was  selected  as  one  of  the  GPS 
positioning  receive!  s  to  acquire  data  for  the  experiment.  One  version  of  the  TI-4100  is 
for  commercial  users  and  the  other  (TI-4100  GEOSTAR)  is  a  tri-agency  version  for  use 
by  the  Defense  Mapping  Agency  (DMA),  National  Oceanic  and  Atmospheric 
Administration  (NOAA),  and  the  United  States  Geological  Survey  (USGS).  TI-4100 
GEOSTAR  software  has  been  obtained  for  use  at  NPS  from  NOAA  and  the  Naval 
Surface  Weapons  Center  (NSWC)  in  Dahlgren,  Virginia,  for  compatibility 
modifications  with  NFS's  IBM  3033  mainframe  computer.  Of  this  software,  three 
programs  from  NSWC,  written  chiefly  by  S.  L.  Mcycrhoff  for  their  CDC  Cybcr-865 
computer,  reformat  GPS  data  and  then  compute  receiver  positions  and  time  by  either 
an  iterative  least  squares  approach  or  a  Kalman  filter. 

The  main  objectives  of  this  thesis  are  to  describe  the  effr  '  modifying  NSWC's 
software  codes  and  algorithms  for  compatibility,  and  establishing  them  on  the 
IBM-3033  mainframe  computer  at  NPS.  In  addition,  discussions  of  the  mathematical 
models,  results  of  test  runs  and  validation  of  the  modified  software  have  been  included. 
This  software  was  acquired  from  NSWC  for  exclusive  use  at  NPS.  Future  users  at 
NPS,  are  encouraged  to  test  and  verify  their  solutions  as  correct  before  accepting  them. 
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II,  BACKGROUND 


A.  GLOBAL  POSITIONING  SYSTEM 

GPS  is  a  universal  satellite  positioning  system  that  is  independent  of  weather, 
time,  or  geographical  position.  It  provides  both  position  and  time  in  an  earth-fixed 
geocentric  cartesian  coordinate  system  as  well  as  velocity  data.  Since  it  is  a  space 
based  radio  navigation  system,  it  can  provide  accurate  information  to  an  unlimited 
number  of  users  globally  on  or  near  the  earth's  surface.  Three-dimensional  position 
data  to  within  a  16-m  spherical  error  probability  (SEP),  velocity  to  within  0.06-0.15  m/s 
and  GPS  system  time  tc  within  25  ns  are  possible  instantaneously  to  properly  equipped 
users  anywhere  within  500  nmi  of  the  earth  (Milliken  and  Zoller,  1980). 

The  GPS  is  comprised  of  three  major  segments,  the  Space,  the  Control,  and  ‘he 

User. 

1.  Space  Segment 

The  GPS  when  fully  operational  will  consist  of  18  active  satellites  (Space 
Vehicles  or  SVs)  in  six  planes  with  near  circular  orbits  at  altitudes  of  10,900  nmi  and 
periods  of  12  hours.  Each  of  the  six  orbital  planes  is  inclined  55  degrees  to  the  equator 
and  is  to  contain  three  SVs  equally  spaced.  The  configuration  of  the  SVs  in  the  orbital 
planes  is  such  that  a  minimum  of  four  SVs  will  be  available  to  the  user  at  all  times, 
anywhere  on  or  near  the  earth's  surface.  Each  individual  SV  transmits  its  broadcast 
ephemeris  and  clock  correction,  plus  an  almanac  which  contains  orbital  parameters  and 
clock  correction  estimates  for  all  other  SVs  in  the  system. 

Each  SV  transnuts  two  radio  frequencies  simultaneously,  LI  and  L2 
respectively,  containing  navigation  information  in  the  navigation  data  message.  The 
LI  frequency  is  centered  on  1575.42  MHz  and  the  L2  frequency  is  centered  on  i227.60 
MIIz.  These  signals  arc  in  a  frequency  range  that  provides  good  all  weather  operation 
and  allows  the  user  to  determine  ionospheric  propagation  delay.  The  LI  signal  is 
modulated  by  a  10.23  MIIz  clock  rate  Precision  (P)  code  and  a  i .023  MIIz 
Course/Acquisition  (C/A)  code  (Spilkcr,  1980).  The  P-codc  is  a  random  sequence 
binary  code  that  provides  the  user  the  capability  for  high-precision  positioning  and  is 
highly  resistant  to  electronic  and  multipath  interference.  The  C/A-code  gives  all  users 
rapid  acquisition  capability  and  acts  as  an  aid  to  gain  access  to  the  P-code.  The 
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navigation  data  message  allows  the  user's  receiver  to  calculate  position,  velocity  and 
GPS  time  by  providing  broadcast  epuemeris,  clock  corrections  and  almanacs  for  all 
SVs.  The  L2  signal  is  modulated  by  the  P-code  only. 

2.  Control  Segment 

The  control  segment  consists  of  a  master  control  station,  a  number  of  monitor 
stations  and  ground  antennae  located  around  the  world.  The  USAF  master  control 
station  is  at  Colorado  Springs,  Co.,  and  monitor  stations  are  located  on  Kwajalcin 
Atoll,  Ascension  Island,  Hawaii  and  Diego  Garcia.  DMA  also  maintains  three  stations 
in  the  United  Kingdom,  Australia  and  Argentina,  and  stations  are  proposed  for 
Ecuador  and  Bahrain.  The  monitor  stations  acquire  ranging  data  for  each  SV  in  a 
passive  mode.  These  data  are  then  relayed  to  the  master  control  station  where  they 
are  processed  to  determine  the  SV's  position  and  signal  data  accuracy  (Wooden,  1985). 
The  master  control  station  takes  the  relayed  data,  computes  errors,  generates  a  new 
navigation  data  message  and  then  uploads  this  message  to  each  SV  via  the  ground 
antennae.  The  ground  antennae  transmit  as  well  as  receive  satellite  contrc. 
information. 

3.  User  Segment 

This  segment  consists  of  stationary,  low,  medium  and  high  dynamic  receivers 
designed  with  specific  requirements  for  individual  users.  The  user  receivers  are 
designed  to  receive  and  process  SV  data  from  four  SVs  either  simultaneously  or 
sequentially.  The  TI-4100  receiver  determines  the  range  from  the  SV  by  measuring  the 
time  delay  of  the  SV's  specific  epoch  and  the  user  generated  codes.  Then  making 
separate  ranging  and  carrier  frequency  Doppler  shift  measurements  to  four  SVs,  the 
receiver  calculates  the  user's  position,  clock  bias  error,  clock  drift  and  velocity.  The 
position  coordinates  are  in  the  World  Geodetic  System  (WGS)  1972,  an  earth-fixed 
geocentric  cartesian  coordinate  system. 

B.  TEXAS  INSTRUMENTS  4100  NAVSTAR  RECEIVER  (TI-4100) 

1.  Characteristics  and  Description 

The  TI-4100  receiver  is  portable  for  field  operations  (it  is  about  the  size  of  a 
small  suitcase)  but  can  be  mounted  in  a  standard  19-inch  rack  if  desired.  A  small 
lightweight  hand-held  Control  Display  Unit,  contains  a  keyboard  and  a  display 
window.  The  receiver  is  modular  in  design  and  allows  for  easy  removal  and 
replacement  of  circuit  boards. 
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The  TI-4100  is  a  single  channel  dual-frequency,  multiplexing  NAVSTAR 
receiver  and  navigation  processor  that  tracks  up  to  four  SVs  (Texas  Instruments, 
1983).  The  dual  frequency  receiver  allows  for  tracking  LI  and  L2  frequencies  for  each 
of  the  four  SVs.  It  also  has  an  optional  dual  cassette  recorder  for  recording  GPS/SV 
data.  If  GPS/SV  data  is  to  be  acquired  in  a  degraded  mode  (i.e.  less  than  four  SVs) 
there  is  the  capability  of  adding  an  external  atomic  clock  time  standard  input  to  the 
TI-4100  receiver.  The  TI-4100  receiver,  Control  Display  Unit,  antenna  and  cassette 
recorder  are  shown  in  Figure  2.1. 

2.  Receiver  Measurements  and  Post  ion  Calculations 

The  basic  measurements  made  by  the  TI-4100  receiver  are  the  arrival  times  of 
the  SV-generated  signals  compared  to  the  receiver  clock,  which  contains  a  stable 
oscillator  and  a  counter  (Texas  Instruments,  1983).  Four  pseudorange  measurements 
plus  the  navigation  data  message  provide  enough  information  to  calculate  a 
three-dimensional  user  position  solution  plus  the  user  clock  correction  to  true  GPS 
time.  The  solutions  are  obtained  using  four  observation  equations  and  four  unknowns. 
Doppler  shifts  due  to  the  relative  velocity  between  the  SV  and  the  user  are  measured. 
Measurements  are  also  made  of  the  phase  differences  for  the  LI  and  L2  signals 
between  the  SV  and  the  user. 

3.  Position  Accuracy 

Error  statistics  for  the  TI-4100  receiver  are  given  in  Table  I.  It  should  be 
noted  for  this  table  that  Operating  Modes  using  time  aiding,  have  a  positive  growth 
rate  error  of  2.0  to  10.0  m/h  (i-sigma)  for  frequency  standards  that  allow  navigation 
calculations  to  estimate  precise  frequency  (Texas  Instruments,  1983).  The  positioning 
accuracies  in  Table  I  were  based  on  the  TI-4100  receiver  observing  pseudoranges, 
Doppler  shift  data,  phase  differences  for  LI  and  L2,  then  computing  position,  time  and 
velocity  solutions.  The  pseudoranges  from  the  TI-4100  receivers  are  the  major 
observables  considered  in  this  thesis. 

C.  PSEUDORANGE  DATA  AND  THEIR  CHARACTERISTICS 

The  basic  operation  of  the  GPS  depends  on  the  user  determining  pseudorange 
and  range  rate  to  a  number  of  GPS  SVs  which  have  precisely  known  cphemcrides 
(Texas  Instruments,  1983).  Pseudorangc  is  the  apparent  range,  which  includes  the 
clock  errors  between  the  SV  and  the  user's  receiver.  The  transit  time  of  the  navigation 
signal  between  each  SV  and  the  user  is  measured  by  the  receiver  and  then  scaled  by  the 
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TABLE  I 

TI-4100  POSITION  ACCURACY 


(From  Texas  Instruments,  TI-4100  Owner's  Manual,  1983) 

1-Sigma  1-Sigma 

Position  Error  Velocity 
Error  (m)  ,  Error 

Operating  Mode  Statistics  P  C/A  (m/s) 


Four  satellites 


Stationary 

0.5  g  acceleration 

2  g  acceleration 

4  g  acceleration 

Differential 

Stationary  differential 

3-D 

3-D 

3-D 

3-D 

3-D 

3-D 

13 

14 

14 

30 

4 

3 

46 
'  47 

47 

55 

12 

8 

N/A 

0. 15 
2.  00 
5.00 
0. 15 
N/A 

Three  satellites 

Degraded  (altitude  aiding) 

Horizontal 

16 

24 

0.15 

Differential  degraded 
(altitude  aiding) 

Horizontal 

5 

6 

0.15 

Degraded  (external  cesium 
frequency  standard  time 
aiding) 

3-D 

16 

47 

0.15 

Differential  Degraded 
(external  cesium  frequency 
standard  time  aiding) 

3-D 

5 

12 

0.15 

Two  satellites 

Degraded  (altitude  and 
time  aiding) 

Horizontal 

16 

24 

0. 15 

Differential  degraded 
(altitude  and  time  aiding) 

Horizontal 

5 

6 

0.15 

speed  of  light  to  compute  a  raw  pseudorange.  Each  raw  pseudorange  must  be 
corrected  for  tropospheric  and  ionospheric  propagation  delays,  SV  clock  offset,  user 
clock  bias,  relativity  corrections,  and  earth  rotation  corrections.  The  pseudorange 
calculation  is  given  by  Equation  (2.1). 

Rj*  =  Rj  +  cAtAi  +  c(Atu  -  Atsi)  (2.1) 

where  i  =  1  to  4  (pertaining  to  four  observed  SV's) 
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r.*  =  pseudorange  to  the  SV 

Rj  =  true  range 

c  =  the  speed  oflight 

Atsi  =  SVj  clock  offset  from  GPS  system  time 

Atu  =  user  clock  offset  from  GPS  system  time 

AtAi  =  propagation  delays  and  other  errors 

The  tropospheric  delay  is  estimated  using  simple  refraction  models  with  standard 
atmospheric  values.  Ionospheric  delay  is  calculated  by  taking  advantage  of  the  dual 
frequency  measurements  of  the  LI  and  L2  signals.  The  SV  clock  error  (difference 
between  SV  clock  and  GPS  system  time)  is  modeled  by  the  Control  Segment  and  these 
correction  coefficients  are  transmitted  in  the  navigation  data  message  from  the  SV  to 
the  user.  User  clock  bias  is  due  to  the  fact  that  the  receiver  clock  is  not  perfectly 
synchronized  to  the  SV  clock.  Relativity  corrections  are  due  to  the  eccentricity  of  the 
SV's  orbit.  Earth  rotation  corrections  are  due  to  the  motion  of  the  earth  during  the 
signal's  propagation  from  SV  to  the  user.  Once  four  raw  pseudoranges  have  been 
corrected  for  systematic  errors,  they  are  used  to  calculate  the  receiver's 
three-dimensional  position  and  time. 

The  precision  of  position  calculated  by  pseudoranges  is  dependent  on  the 
geometry  of  the  SV's  with  respect  to  the  user.  This  geometry  is  continually  changing, 
even  if  the  user's  receiver  is  not  in  motion  relative  to  the  earth.  If  the  intersection  of 
the  pscudoranges  is  orthogonal,  errors  in  the  pseudorange  measurement  are  minimized 
whereas  if  the  angle  of  intersection  is  small,  errors  are  large  because  of  the  geometric 
relationship. 

Geometric  Dilution  of  Precision  (GDOP)  is  a  composite  "indicator"  expressing 
the  influence  of  the  SV's  geometry  on  the  accuracy  of  the  user's  clock  offset  and 
position  determinations  (Milliken  and  Zoller,  1980).  Other  terms  to  describe  dilution 
of  precision  in  position  and  time  are:  Position  Dilution  of  Precision  (PDOP)  which 
yields  dilution  of  precision  in  three-dimensions;  Horizontal  Dilution  of  Precision 
(HDOP)  which  gives  dilution  of  precision  in  the  two  horizontal  dimensions;  Vertical 
Dilution  of  Precision  (VDOP)  which  yields  the  dilution  of  precision  in  the  vertical 
dimension;  and  finally  Time  Dilution  of  Precision  (TDOP)  is  the  dilution  of  precision  in 
the  user's  time  bias  (Milliken  and  Zoller,  1980).  GDOP,  which  is  dimensionless,  is 
calculated  by  Equation  (2.2),  (Milliken  and  Zoller,  1980), 
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(2.2) 


GDOP  =  [(PDOP)2  +  (TDOP)2]1/2 

The  best  geometry  for  a  four  SV  configuration  is  the  combination  which  results  in  the 
lowest  GDOP. 

Once  pseudorange,  Doppler  shift  and  phase  data  have  been  acquired,  and  the 
receiver's  field  solutions  calculated,  post-processing  of  these  data  may  yield  improved 
position  and  time  solutions.  To  make  TI-4100  data  accessible  to  other  users  with 
different  computer  processing  systems,  it  has  become  necessary  to  format  SV  data  in  a 
manner  that  can  be  exchanged  among  independent  users. 

D.  STANDARDIZED  EXCHANGE  FORMAT  FOR  NAVSTAR  GPS 

GEODETIC  DATA 

This  Standardized  Exchange  Format  (SEF)  was  proposed  by  V.  Dan  Scott  and  J. 
Gary  Peters  (1983),  and  was  adopted  to  make  possible  the  universal  exchange  of 
NAVSTAR  GPS  data.  To  do  this,  serious  compromises  were  made  in  regard  to 
optimum  magnetic  tape  utilization,  computer  execution  efficiency,  and  computer 
program  size  (Scott  and  Peters,  1983).  These  compromises  were  required  to  allow  for 
maximum  tape  compatibilty  between  different  computer  systems. 

1.  The  Basic  Format 

The  building  block  for  each  standardized  data  exchange  tape  is  the  record, 
which  is  80  characters.  Under  this  format,  an  unlabeled  9 -track,  ASCII-coded  tape 
with  an  80-character  blocking  factor  is  used.  The  header  (or  the  first)  file  on  a  tape 
consists  of  comment  records  and  is  not  considered  part  of  the  data  set  but  is  used  to 
help  the  user  to  identify  the  contents  of  the  tape.  When  the  last  data  file  is  written,  on 
the  tape  a  special  end  of  information  file  is  written  (Scott  and  Peters,  1983). 

On  the  tape,  records  fall  into  two  categories:  Control  Information  or  Data 
Information  records.  Control  Information  records  are  used  to  identify  which 
FORTRAN  Format  statement  to  use  when  reading  data  records.  Data  Information 
records  consist  of  one  of  eight  data  items  that  are  uniquely  identified  by  a  four  digit 
data  item  number.  The  eight  data  items  under  the  SEF  are  listed  below: 

(1)  Constant  equipment  and  software  specific  data 

(2)  Campaign  and  station  specific  data 

(3)  Calibration  data 

(4)  Time  tag  update  data 

(5)  Space  Vehicle  specific  data 


16 


(6)  Weather  data 

(7)  Real  time  solution  data 

(8)  Measurement  data 

The  goal  of  this  format  is  to  be  independent  of  hardware  and  software  design.  In 
addition  to  data  formats,  the  reference  Coordinate  System  for  the  data  must  be  defined. 
2.  Coordinate  System 

The  WGS  72  an  earth-fixed  geocentric  cartesian  coordinate  system  is  defined 
as  follows  (Seppelin,  1974): 

X-Axis  =  Intersection  of  the  WGS  72  reference  meridian  plane 
and  the  plane  of  the  Equator. 

Y-Axis  =  Measured  in  the  plane  of  the  Equator  90  degrees 
east  of  the  X-Axis. 

Z-Axis  =  Parallel  to  the  direction  of  the  conventional 

international  origin  (CIO)  for  polar  motion. 

The  reference  ellipsoid  for  WGS  72  is  defined  by  the  parameters  semi-major  axis,  a  = 
6,378,135  m,  and  flattening,  f  =  1/298.26  (Seppelin,  1974). 
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III.  DATA  PROCESSING 


The  pre-processing  of  GPS  data  is  carried  out  to  obtain  a  format  suitable  for 
further  manipulation.  This  is  necessary  when  geodetic  positioning  of  the  receiver 
(either  static  or  dynamic)  is  desired,  through  post-processing  the  SV  data.  Figure  3.1 
portrays  the  flowchart  for  processing  the  TI-4100  data. 

A.  PRE-PROCESSING 

1.  Manipulation  of  Cassette  Recorded  Data 

The  first  step  in  processing  GPS  data  recorded  on  cassettes,  is  to  dump  the 
raw  data  from  the  cassettes  to  either  a  9-track,  magnetic  tape  or  disk.  The  Applied 
Research  Laboratories  (ARL)  at  Austin,  Texas  have  developed  software  to  accomplish 
this  fundamental  task.  The  software  at  ARL  can  dump  and  process  the  raw  data 
recorded  on  cassettes,  and  format  the  GPS  data  onto  a  9-track  magnetic  tape  in  the 
SEF  (Chapter  II. D.).  This  tape  can  then  be  transferred  to  disk  for  convenience  in 
further  processing  with  a  mainframe  computer. 

2.  Alignment  of  Common  Time-Tagged  Data 

The  alignment  of  common  time-tagged  data  items  is  accomplished  through  use 
of  the  CON9TR  program. 

a.  Program  Description 

Program  CON9TR,  the  first  of  three  programs  acquired  from  NSWC 
(Chapter  I.),  reads  the  SEF  data  and  formats  it  into  another  NSWC  specified  file.  This 
file  then  becomes  the  input  for  the  two  remaining  programs,  SOLTOM  and  KALMN2 
and  has  all  common  time-tagged  data  on  the  same  record.  The  only  calculations 
performed  at  this  stage  are  variances  on  the  pseudorangc  data  and  the  delta  range 
(Doppler)  data  which  are  measures  of  the  signal  to  noise  ratio. 

The  program  also  looks  at  the  "quality  vectors"  (Texas  Instruments, 
1982)  of  the  data  on  the  9-track  SEF  tape  and  returns  an  integer  value  representative 
of  the  data  quality  at  that  time  line.  Four  integer  values  are  possible: 

0  =  good  data 

1  =  good  data  but  this  is  the  first  Doppler  count 

14  =  bad  data  at  that  point  but  there  was  not  a  break 

in  the  Doppler  data 

15  =  bad  data 
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Figure  3.1  Flowchart  For  Processing  TI-4100  GBOSTAR  GPS  data. 
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The  C0N9TR  program  requires  a  user  input  file  (Figure  3.1)  in  addition 
to  the  GPS  data  in  the  Standardized  Exchange  Format  during  program  execution. 
CON9TR  consists  of  the  main  program  and  22  subroutines.  Appendix  A  gives  a  brief 
description  of  each  subroutine's  function  and  nesting. 

b.  Input/Output 

The  user  input  file  (Figure  3.1)  for  the  CON9TR  program  consists  of  the 

following: 

(1)  The  number  of  trackers,  usually  four  but  may  be  less  than  four 

(2)  GPS  time  of  the  first  data  point  to  be  put  on  the  new  file 

(3)  Time  to  end  data  on  the  new  file,  if  the  end  of  the  file 

is  found  before  the  time  to  end  data,  then  data  will  end 
at  the  end  of  file 

(4)  Range  variance  bias  factor  [km^] 

(5)  Doppler  variance  bias  factor  [km^] 

(6)  X-coordinate  of  the  initial  receiver  position  that 
will  be  put  on  the  new  files 

(7)  Y-coordinate  of  the  initial  receiver  position  that 

will  be  put  on  the  new  files 

(8)  Z-coordinate  of  the  initial  receiver  postion  that 

will  br  put  on  the  new  files 

Appendix  B  gives  the  user  input  guide  to  CON9TR  program  variables  and  format. 

The  first  output  file  from  CON9TR  (Figure  3.1)  writes  all  common  time 
satellite  tracking  data  into  the  same  record  block,  and  Appendix  C  gives  the  output 
format.  The  second  output  file  writes  out  messages  to  aid  the  user  in  data  quality 
analysis,  flagging  one  of  four  quality  vector  conditions  that  exist  in  the  data. 

c.  Compution  of  Variances 

The  variances  for  the  pseudorange  and  Doppler  data  in  the  CON9TR 
program  are  computed  by  Equations  (3.1)  and  (3.2),  respectively  (Texas  Instruments, 
1982). 

«2codc  -  KbL  *  N)/(IO(CVIObl  *  (3.1) 

(0.5  +  BpD/(IO<Ci/10))]  *  Kj  +  K2 

where  the  asterick  ,*,  is  used  to  signify  multiplication 
Bj^  =  code  loop  noise  bandwidth  Hz 
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Bpo  “  code  loop  predetection  bandwidth  Mz 
Cj  =  Lj  carrier  power  to  noise  spectral  density  ratio 
dB-I-Iz 

N  =  tracking  mode,  1,  2,  3,  4  based  on  SV  mode 
9  9 

Kj  =  (chips)'1  to  [nr]  conversion: 

(29.305)2  for  P-code 
(293.05)2  for  C/A-code 
K2  =  variance  bias  factor  [nr] 

cr2pLL  =  PL  *  N)/(10(Ci/10))]  *  {1  +  [0.5  *  (3.2) 

BpD/(10(Ci/10))]}  "  K3  4-  K4 

where, 

PLL  =  carrier  phase-locked  loop 

Bj^  =  carrier  loop  bandwidth  Hz 

BpD  =  carrier  predetection  bandwidth  Hz 
9  9 

K-j  =  [rad“]  to  [nr  ]  conversion 
=  (c/2Jtf)2 

c  -  speed  of  light 
f  =  L  band  carrier  frequency 
Lj  =  1.57524  x  109  Hz 
L2  =  1.2276  x  109  Hz 

=  variance  bias  factor  [nr] 

d.  Program  Modifications  at  NPS 

The  original  code  from  NSWC  was  written  in  FORTRAN  IV.  This 
code  was  modified  at  NPS  to  be  compatible  with  the  FORTRAN  77  compilers 
(FORTVS  and  WF77)  on  the  IBM-3033  mainframe  computer.  Changes  to  the  original 
code  included  the  following, 

•  Changes  in  the  formatted  write  statements 

•  Changes  in  variable  alignments  in  the  common  statements 

(i.e.  making  sure  double  precision  variables  arc  aligned 
ahead  of  single  precision  variables) 

•  Changes  in  dimension  size  of  the  data  arrays 

•  Changes  in  defining  variables 
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•  Changes  in  the  placement  of  data  intialization 

statements 

•  Addition  of  the  Block  Data  subroutine  to  contain  the 

moved  data  intialization  statements  and  accompanying 
common  statements. 

Executive  files,  written  at  NPS  to  handle  t.ne  input  and  output  files  for 
the  CON9TR  program,  included  the  auto-aouble  precision  commands  for  all  real 
variables  using  the  FORTVS  compiler.  This  modification,  automatically  doubles  the 
declared  precision  in  the  code  and  was  required  due  to  ihe  use  of  60  bit  words  by  the 
CDC  Cyber-865  computer  and  the  32  bit  words  by  the  IBM-3033  mainframe  computer 
at  NPS. 

B.  DETERMINATION  OF  PSEUDORANGE  POSITION  SOLUTION 

User  receiver  position  solutions  whether  fixed  (static)  or  moving  (dynamic)  may 
be  calculated  by  post-processing  the  TI-4100  GEOSTAR  GPS  data. 

1.  Static  Positioning 

The  second  program,  SOLTOM  from  NSWC  (Chapter  I.,  computes  a  position 
solution  for  a  static  receiver  by  using  an  iterative  batch  least  squares  calculation. 

a.  Program  Description 

Program  SOLTOM  computes  receiver  position  and  time,  solving  for 
initial  corrections  to  estimates  of  parameters  X,  Y,  Z,  time  biases  and  scaling  factor  for 
a  tropospheric  correction.  Types  of  solutio'  :  which  can  be  selected  through  the  user 
input  file  are  described  in  Chapter  III.B.i.c..  The  program  consists  of  the  main 
program  and  52  subroutines.  Appendix  D  gives  a  brief  description  of  each  subroutine's 
function  and  nesting. 

b.  Math  Models 

The  math  models  for  the  batch  least  squares  solution  in  program 
SOLTOM  solve  for  corrections  to  the  initial  estimates  of  the  parameters  desired.  The 
observation  equation  is  given  by  Equation  (3.3): 

R°i  =  l(Xsv-  Xp)2  +  (YSV-  Yp)2  +  (Zsv-  Zp)2]1/2  +  CB  +  CD  +  CA  +  TR  (3.3) 
where, 

R°;  =  i^1  range  observed  from  SV  of  receiver 
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Xsy  =  earth-fixed  geocentric  X-coordinate  of  SV 
Ygv  =  earth-fixed  geocentric  Y-coordinate  of  SV  ■ 

Z$v  =  earth-fixed  geocentric  Z-coordinatc  of  SV 
X  !=  earth-fixed  geocentric  X-coordinate  of  receiver 

F 

Y  =  earth-fixed  geocentric  Y-coordinate  of  receiver 

F 

Zp  •■=  earth-fixed  geocentric  Z-coordinate  of  receiver 

CB  =  Clock  bias  term,  SV  or  receiver  depending  on  user 
input  selection 

CD  =  Clock  drift  term,  SV  or  receiver  depending  on  user 
input  selection 

CA  =  Clock  aging  term  (frequency  drift),  SV  or  receiver 
depending  on  user  input  selection 

TR  =  Tropospheric  refraction  correction. 

The  least  squares  solution  of  Equation  (3.3)  in  matrix  form,  can  be  written  as: 

BP  +  E  =  0  (3.4) 

where, 

B  -  Matrix  containing  the  partial  derivatives  of  the  obser  ation 
equation  with  respect  to  the  parameters  being  solved  for 

P  =  Vector  containing  the  differences  of  the  parameter  adjusted 
(Pa)  and  the  parameter  estimate  (Pe),  (Pa  -  Pe) 

E  -  Vector  containing  the  partial  derivatives  of  the  observation 
equations  multiplied  by  the  weight  matrix,  multiplied  by  tuc 
data  residuals 

The  matrix.  B  is  obtained  from  Equation  (3.5), 

B  =  ATwA  (3.5) 


23 


where, 

A  =  Partial  derivatives  of  data  equation  with  respect  to  a 
parameter  j 

AT  =  Transpose  of  the  A-matrix 

w  =  Weight  matrix  with  on  ^ie  diagonal,  zeros 

off  diagonal,  where  k  =  1  to  number  of  data  points,  this 
is  the  variance  of  noise  level  on  a  "k"  range  measurement 

or  writing  Equation  3.5  in  summation  form. 

By  =  Enk=  1  /  api  *  aR°k  /  apj  *  *k2  <3-«> 

where, 


n  =  number  of  data  points  used  in  the  least  squares  solution 

R^o  =  k^  data  equation 

Pj  =  i^1  parameter  being  solved  for 

Pj  =  j1*1  parameter  being  solved  for 

<r^2  =  k^  variance  vaiue 

The  A-matrix  or  design  matrix  is  defined  by  Equation  (3.7): 

Ay  -  3R°i  /  aPj  (3-7) 

where, 

r°.  _  jth  (jata  p0jn(;  (3.3) 

j» 

Pj  =  j'n  parameter  being  solved  for 

The  inverse  of  B  is  the  variance-covariance  matrix,  COV,  with  the  variances  <sc  of  the 
parameter  corrections  AP  that  are  being  solved  for,  on  the  diagonal  given  by  Equation 
(3.8), 
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B'1  =  COV 


(3.8) 


where,  COV  is  a  6x6  matrix  if  six  parameters  are  being  solved,  or  7  x  7  if  the  optional 
7^  parameter  (CA)  is  to  be  determined.  The  off  diagonal  elements  are  the 
cross-correlations  of  the  paired  parameters  for  that  time  line.  The  vector  E  in 
summation  form  is  shown  in  Equation  (3.9), 

Ei  “  IV  1  <3Rk°/a  ‘V  *  V  *  Rres 

where, 

R^o  =  k^1  data  equation 

n  =  number  of  data  points  used  in  least  squares 
calculation 

=  i^1  parameter  being  solved  for 
<Tj^  =:  k^1  weighted  variance  value 
Rres  =  range  residuals 
or  writing  Equation  (3.9)  in  matrix  form, 

E  -  RreswA  (3.10) 

where, 

AP  =  B'*E  (3.11) 

and 

AP  =  corrections  to  the  initial  estimates  of  unknown  parameters  being  solved  for  (X, 
Y,  Z,  clock  biases,  and  tropospheric  refraction  correction). 

The  following  Equations  (3.12)  to  (3.35)  describe  the  various  corrections,  Time  Epoch, 
clock,  Ionospheric  refraction,  Tropospheric  refraction,  earth  rotation,  and  relativistic 
effects  which  are  required  to  correct  pscudorange  and  delta  range  (Doppler)  data 
(Mcycrhoff,  1985). 
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Time  Epoch  Correction 


Tsv-tfT0c  (312) 

where, 

Tsv  =  time  from  epoch 

tt  =  the  GPS  time-tag  for  range  and  delta  range  (Doppler) 
data 

Toc  =  time  of  reference  for  clock  corrections 

To  correct  the  1'  for  the  end  of  the  week  cross  over,  correct  T$v  as  follows: 

If  T$v  <  302400  [s]  then,  Tsv  =  Tsy  +  604800  Is] 

If  Tsv  >  302400  [s]  then,  Tsv  =  T$v  -  604800  [s] 

Clock  Correction 

ccl(tt)  =  [A0  +  (A,  *  Tsv)  +  (Aj  *  Tsv  *  Tsv)]  *  c  (3.13) 


ACcl(U)  =  [A0  +  (A,  *  {(2  *  Tsv)  -  td))]  *  td  *  c  (3.14) 

where, 

Cc|(tt)  =  clock  correction  for  pseudorange  data 
at  GPS  time-tag  in  [km] 

ACcj(tt)  =  clock  correction  for  delta  range 

(Doppler)  data  at  GPS  time-tag  in  [km] 

A0  =  clock  correction  bias  term  [s] 

Aj  =  clock  correction  drift  term  [s/s] 

A2  =  clock  correction  aging  term 
[s/(s)2| 

T$v  =  time  from  epoch 
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c  =  speed  of  light  299792.458  [km/s] 

td  =  the  Doppler  count  interval  for  delta  range 
(Doppler)  data 

Ionospheric  Refraction  Correction 

Ci0(tt)  =  CPj(tt)  -  CP2(tt)/(Qj/Q2)2  ■  1.0  (3.15) 


DCi0(tt)  =  [(DQPt  , (tt)/Q , )  -  (DOPf  2(tt)/Q2)]  *  c  (3.16) 

[(Qi/Q2)2  - 1.0)  vos 

where, 

Ci0(tt)  =  pseudorange  ionospheric  refraction 
delay  in  [km] 

tt  =  GPS  time-tag  for  pseudorange  and  delta  range 
(Doppler) 

CPj  =  LI  pseudorange  in  [km] 

CP2  =  L2  pseudorange  in  [km] 

DOPli  =  LI  Doppler  count 

DOPl2  =  L2  Doppler  count 

DCj0  =  delta  range  (Doppler)  frequency  shift  due 
to  ionosphere 

Qj  =  LI  frequency  multiplier  (154.0) 

Q2  =  L2  frequency  multiplier  ( 1 20.0) 

VQS  =  The  nominal  satellite  oscillator  frequency,  10.23  MHz 

In  Equation  (3.15)  the  correction,  Cj0(tt),  is  used  to  correct  the  pseudorange  that  is 
time-tagged,  tt.  In  Equation  (3.16)  the  correction,  DCj0(tt),  is  applied  to  the  delta 
range  (Doppler)  data  to  correct  it  for  time-tag,  tt. 
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Tropospheric  Refraction  Correction 

The  Tropospheric  refraction  computation  can  be  divided  into  two  parts, 
weather  dependent  and  receiver  geometry  ..alculations  (Meyerhoff,  1985). 

Weather  Dependent  Correction 

First,  convert  temperature  to  degrees  Kelvin  and  convert  humidity  to  a 

fraction: 

Tk  =  Tc  + 273.0  (3.17) 


RH  =  Hd/ 100.0  (3.18) 

where, 

Tc  =  air  temperature  in  degrees  centigrade 
Hd  =  the  relative  humidity  in  percent. 

Second,  compute  the  surface  partial  pressure,  E0,  from  temperature  and  relative 
humidity. 

E0  =  RH  *  35.65  *  10  exp[7.617-(2285.0/Tk)]  (3.19) 

Third,  the  components  of  the  zenith  value  of  tropospheric  delay,  Zdry  and  Zwet  are 
calculated  from,  EQ,  and  atmospheric  pressure. 

Zdry  =  2.276  *  Pb  *  0.01  (3.20) 

Zwct  =  lm  *  E0cxp(1.23)/Tk)  +  (1.705  *  106  (3.21) 

*  Ah  *  E0oxp(I.46/Tk3) 

where, 

exp  =  exponent 

Pb  =  the  atmospheric  pressure  in  [mb] 

Ak  =  temperature  lapse  rate  set  at  0.006  °C/m 
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exp  =  exponent  It  should  be  noted  that  Zj^  and  Zwet  need  only  be  recomputed, 
when  the  weather  data  changes. 

Receiver  Geometry  Dependent  Correction 

First,  to  find  the  elevation  angle,  EV(t),  of  the  satellite  at  time,  t,  take 
the  inverse  sine  of  the  dot  product  of  the  unit  vector  from  the  receiver  to  the  satelf  e 
and  the  unit  position  vector  of  the  receiver. 

EV(t)  =  sin-1  RS*  ~R  (3.22) 

MGrs  *  MG^ 

where, 

RS  =  vector  from  the  receiver  to  the  satellite 
MGrs  =  vector's  magnitude 
R  =  position  vector  of  the  receiver 
MGr  =  magnitude  of  R 


MGrs  =  [(X(t)  -  Xr)2  +  (Y(t)  -  Yf)2  +  (Z(t)  -  Zj.)2]1/2  (3.23) 


M  Gr  =  [Xr2  +  Yr2  +  Zj.2]1/2 


(3.24) 


where, 

Xr,  Yf,  Zr  =  the  estimate  of  the  receiver's 
position  in  an  earth-fixed 
geocentric  coordinate  system 

rewriting  Equation  (3.23),  EV(t)  can  be  written  as, 

SEV(t)  =  [{(X(t)  -  Xr)  *  Xr)  (3.25) 

+  {(Y(t)  -  Yr)  *  Yr)} 

+  {(Z(t)  -  Zr)  *  Zr})/(MGrs  *  MGr) 

where, 

X(t),  Y(t),  Z(t)  =  the  position  coordinates  in  a  earth- 
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fixed  geocentric  coordinate  system 


or 

EV(t)  =  sin1  [SEV(t)]  (3.26) 

Second,  the  multiplication  factors  for  dry  and  wet  as  a  function  of  elevation  are 
calculated  from  the  elevation  angle,  EV,at  GPS  time- tag  (tt). 

Fdry(t)  =  1  /  [sin  EV(t)J  +  {0.00 143/[( tan  EV(t))+  0.0445]}  (3.27) 

Fwet(t)  =  1  /  [sin  EV(t)]  +  {0.00035/[(tan  EV(t))+  0.17]}  (3.28) 

Third,  Fdry(t  -  td)  and  Fwct(t  -  td)  are  calculated  in  a  similar  manner  as  Equations 
(3.27  and  3.28). 

Total  Tropospheric  Refraction  Correction 

TR(tt)  -  lZdry  *  Fdry(t)]  +  [Zwet  ♦  Pwc,(t)J  (3.29) 

DTR(tt)  -  <Zdry  *  lFdry(t)  -  Fdry(t  -  td)])  +  (Zwc,  *  |Fwet(t)  -  (3.30) 

FwetO ' “>® 

where, 

TR(tt)  and  DTR(tt)  arc  range  corrections  are  in  [km]. 

Earth  Rotation  Correction 

An  earth  rotation  correction  is  required  due  to  the  rotation  of  the  earth 
during  the  signal  propagation. 

Ccr  =  [(Y(t)  *  Xr)  -  (X(t)  *  Yr)]  *  RRE/c  (3.31) 
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DCcr  =  (({Y(t)  -  Y(t  -  td)}  *  Xr)  -  ((X(t)-X(t  -  td)}  *  Yr)]  *  RRE/c 


(3.32) 


where, 

X(t),  Y(t),  Z(t)  =  position  coordinates  in  an  earth- 
fixed  geocentric  cartesian 
coordinate  system  at  pseudorange 
transmit  time  "t"  used  with  data 
time  tagged  "tt" 

Xf,  Yr  =  estimate  of  the  X  and  Y  coordinate 
in  an  earth-fixed  geocentric 
cartesian  coordinate  system 

RRE  =  Rotational  rate  of  the  Earth  corrected  for  motion 
of  the  equinox 
(0.00007292115855  [rads/s] 

td  =  Doppler  count  interval  for  delta  range 
(Doppler)  data. 

Relativistic  Correction 

The  correction  for  the  relativistic  effects  is  due  to  the  eccentricity  of  the 
satellite's  orbit. 

Rel(tt)  =  (2  *  [GM]!/2/c)  *  e  *  A  *  sin  [Ek(t)]  (3.33) 

at  the  beginning  of  the  Doppler  count, 

DRel(tt)  -  (2  *  [GM]1^)  *  e  *  A  *  sin  [Ek(t  -  td)]  (3.34) 

which  yields, 

DRC(tt)  =  RC(t)  -  RC(t  -td)  (3.35) 

where, 

tt  =  GPS  time-tag  for  both  pseudorange  and  delta  range 
(Doppler)  data 
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td  =  Doppler  count  interval  for  delta  range  (Doppler) 
data 

t  =  GPS  time-tag  at  transmission  time  for  data  time- 
tagged  tt 

e  =  eccentricity  of  the  satellite  orbit 

A  =  the  square  root  of  the  semi-major  axis  of  the 
satellite's  orbit 

Ek(t)  =  eccentric  anomaly  for  the  satellite  at  pseudorange 
transmission  time  t 

Ek(t  -  td)=  eccentric  anomaly  for  the  satellite  at  the 
beginning  of  the  Doppler  count  interval 

GM  =  WGS  72  value  of  the  Earth's  gravitational  constant  including 
atmosphere  398600.8  knvfys^ 

RC  =  computed  range 

these  arc  the  defined  constants  and  variables. 

c.  Input/Output 

There  are  two  input  files  for  the  Soltom  program  (Figure  3.1).  The  first 
file  is  the  user  input  file,  where  the  user  sets  which  types  of  solutions  will  be  calculated 
and  which  data  solutions  will  be  printed  out.  The  user  input  guide  for  this  file  is 
described  in  Appendix  E.  The  second  input  file  for  SOLTOM  is  an  output  file  from 
the  CON9TR  program  containing  records  that  are  aligned  with  common  time-tagged 
GPS  data  items. 

There  arc  also  two  output  files  for  SOLTOM  (Figure  3.1).  The  first  file 
is  the  updated  receiver  position,  time  and  other  print  options  such  as  solution  sigmas 
and  clock  bias  terms,  (as  determined  by  the  user  input  file).  The  second  output  file  is  a 
data  corrections  residual  file,  which  can  be  used  in  plotting  routines  and  for  input  to 
other  software  that  handle  multiple  station  positioning  (Mcyerhoff,  1986). 

d.  Program  Modifications  at  NPS 

The  source  code  from  NSWC  for  the  SOLTOM  program  was  written  in  FORTRAN 
IV  which  has  been  modified  to  run  at  NPS  on  the  IBM-3033  mainframe  computer. 
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Similar  changes  to  the  source  code  in  the  CON9TR  program  were  also  made  to  the 
SOLTOM  program  including  the  executive  files  (Chapter  III.A.2.).  The  original  code 
was  written  to  be  able  to  handle  the  precise  ephemeris  data  as  well  as  the  onboard 
broadcast  ephemeris  data  in  the  navigation  data  message  from  the  SV.  The  test  data 
for  SOLTOM  was  generated  at  NSWC  using  the  Broadcast  ephemeris  only,  and  as 
such,  the  modified  program  version  at  NFS  has  the  Precise  Ephemeris  subroutines  and 
corresponding,  call,  statements  commented  out. 

e.  Solution  Techniques  Through  Least  Squares 

There  are  ten  types  of  solutions  possible  with  the  SOLTOM  program 
which  are  dependent  on  the  user  input  file  (Meyerhoff,  1985).  Appendix  E  gives  the 
details  for  the  type  of  solution  desired.  Of  these  ten  types,  nine  are  batch  least  squares 
solutions  and  one  solution  type  five  is  a  Newton-Raphson  solution  (a  geometrical 
solution,  not  a  least  squares  fit).  The  type-five  solution  can  be  used  for  cither  static  or 
dynamic  positioning  performing  an  independent  solution  at  each  time  line.  The  nine 
batch  least  squares  solution  types  fall  into  three  data  classes,  Range  (pscudorange), 
Range  Difference  (delta  range  Doppler),  and  Bias  Range. 

The  reader  is  WARNED ,  that  the  Bias  Range  (solution  type  ten)  has  still 
not  been  tested  by  NSWC,  and  thus  it  should  not  be  used  at  NPS.  Solution  types  one, 
four,  six,  and  nine  are  Range  (pseudorange)  solutions.  Solution  types  two  and  seven 
are  Range  Difference  (Doppler)  solutions.  Finally,  solution  types  three  and  eight  are 
combination  solutions  of  the  Range  (pseudorange)  and  Range  Difference  (Doppler) 
data. 

The  SOLTOM  program  has  the  capability  to  run  more  than  one 
iteration  of  batch  least  squares,  if  needed  (e.g.  for  bad  initial  receiver  position  or  time 
bias  problems).  The  user  has  control  via  the  user  input  file  to  set  how  often  the  least 
squares  solution  will  be  calculated  and  printed  out.  When  a  solution  is  to  be 
computed,  the  program  takes  the  B-matrix  containing  the  partial  derivatives  of  data 
equations  with  respect  to  the  parameters  being  solved  for,  and  multiplies  by  the 
variances  of  the  range  measurements  (Equations  3.5  and  3.6),  inverts  the  B-matrix, 
which  is  the  variance-covariance  matrix  with  variances  of  the  corrections  AP  on  the 
diagonal,  and  computes  the  corrections  to  the  initial  estimates  of  the  parameters  being 
sought,  AP  (Equation  3.11).  Then  Equation  (3.36)  gives  the  updated  station  position 
coordinates,  clock  bias,  drift,  aging,  and  Tropospheric  refraction  correction: 

X  =  X  +  APX  (3.36) 
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Y  =  Y  +  APy 

z  =  z  +  ap2 

~  A0  ^  ^Ao 

Ajt  =  Ajt  +  APAjt 

A2t2  =  A2t2  ■+  APA^ 

TR  =  TR  +  APtr 
where, 

X,  Y,  Z  =  updated  earth-fixed  geocentric  cartesian 
coordinates  through  batch  least  squares 

Aq  =  clock  correction  bias  term  in  [s] 

Ai  =  clock  correction  drift  term  in  [s/s] 

A2  =  clock  correction  aging  term  in 
[s/(s)2] 

(optional  7th  parameter  to  solve  for) 

TR  =  tropospheric  refraction  correction 

The  new  solution  takes  all  previous  solutions  and  averages  up  to  the  current  time  line, 
before  producing  the  iterated  solution. 

2.  Dynamic  Positioning 

The  third  program,  KALMN2  from  NSWC  computes  a  position  solution  for  a 
moving  (dynamic)  receiver,  using  the  UDU'  factorization  of  an  eight-state  extended 
Kalman  filter  (Chapter  I.). 

a.  Program  Description 

Program  KALMN2  uses  the  first  four-states  for  receiver  position  (X,  Y, 
Z)  and  receiver  clock  bias,  and  the  other  four-states  for  receiver  velocity  along  three 
axes  (X,  Y,  Z)  and  receiver  clock  drift.  By  using  the  Doppler  phase  (change  in  range) 
measurements  to  smooth  the  pscudoranges,  signal  multipath  and  measurement  noise 
effects  can  be  suppressed  significantly  to  produce  a  viable  solution  (Mcycrhoff  and 
Evans,  1986). 
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Note:  The  four-states  of  the  velocity  and  dock  drift  are  still  under 
testing  at  NSWC. 

The  KALMN2  program  consists  of  the  main  program  and  43 
subroutines.  Appendix  F,  gives  a  brief  description  of  each  subroutine's  function  and 
nesting. 

b.  Math  Models 

The  math  model  for  the  UD.U'  factorization  of  the  inverse  of  the 
B-matrix  (which  contains  the  partial  derivatives  of  the  observation  equations  with 
respect  to  the  parameters  being  solved  for)  used  in  program  KALMN2  is  described  in 
Bierman  (1977)  and  MeyerhofT  (1986). 

The  TI-4100  receiver  provides  two  pseudorangc  and  two  Doppler 
measurements  (via  LI  and  L2  frequencies),  for  as  many  as  four  satellites 
simultaneously.  The  following  Equations  (3.37  to  3.42)  describe  the  various  corrections 
and  conversions;  receiver  clock  bias,  continuous  count  (Doppler)  to  interval  (Doppler) 
data,  and  Ionospheric  refraction,  which  are  required  for  the  Doppler  and  the 
pseudorange  data  (for  the  Ionospheric  correction)  (MeyerhofT,  1985): 

Receiver  Clock  Dias  Correction 

tt  =  ttr  -  RCB  (3.37) 

where, 

ttr  =  receiver  clocks  GPS  time-tag  for  range  and  delta 
range  (Doppler)  data 

tt  =  GPS  time-tag  corrected  for  receiver  clock  bias 

RCB  =  estimate  of  the  receiver  clock  bias  which  is 

initially  set  to  zero  and  an  improved  estimate 
of  "RCB"  is  made  after  each  time  line 

The  time  delay  data  arc  multiplied  by  the  speed  of  light  "c"  and 
converted  to  a  pseudorange  in  [km]. 

Continuous  Count  to  Interval  ( Doppler )  Data  Conversion 

DOPL1(tt)j  =  DCj(tt)j  -  DCj(tt  -  td)j  -  (td  *  OFFSET!)  (3.38) 
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DOPL2(tt)j  =  DC2(tt)j  -  DC2(tt  -  td)-,  -  (td  *  0FFSET2) 


(3.39) 


where, 

td  =  Doppler  count  interval  for  delta  range  (Doppler)  data 

DCj(tt  -  td) ,  DC2(tt  -  td)=  LI  and  L2  continuous 

Doppler  count  data  at  the  delta  range  interval 

DCj(tt) ,  DC2(tt)  =  LI  and  L2  continuous  count  Doppler 
data  at  the  end  of  the  delta  range  interval 
OFFSET j  ,  OFFSET2  =  Doppler  offsets  used  by  the  receiver 
for  LI  and  L2,  -6000  and  +7600  respectively 

DOP^i(tt) ,  DOPj^2(tt)  =  LI  and  L2  interval  Doppler  data 
for  time-tag  "tt" 

i  =  satellite  from  which  data  originates 

The  delta  range  (Doppler)  data  is  obtained  by  converting  the  LI  interval  Doppler  data 
for  each  satellite. 

RD(tt):  =  (-DOPLU(tt)/V0S  *  Qj)  *  c  (3.40) 

where, 

VQS  =  nominal  satellite  oscillator  frequency,  10.23  MHz 
Qj  =  LI  frequency  multiplier  (154) 

RD(tt)  =  delta  range  (Doppler)  for  one  satellite  in  [km] 

Ionospheric  Refraction  Correction 

RCI(tt)j  -  CPI(tt)j  +  CIO^  (3.41) 

RDCI(tt)j  =  RD(tt)j  -  DCIO(tt)j  (3.42) 

where, 
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CPI  =  uncorrected  pseudorange  measurement  in  [km] 

RD(tt)  =  uncorrccted  delta  (Doppler)  range  measurement 

CIO(tt)  =  ionospheric  refraction  correction  in  [km] 

for  pseudorange  data  (refer  to  Equation  3.15) 

DCIO(tt)  =  ionospheric  refraction  corrections  in  [km] 

for  delta  (Doppler)  range  data  (refer  to  Equation 
3.16). 

Smoothing  For  Noise 

The  pseudorange  data  are  contaminated  due  to  signal  multipath  and 
other  noise  sources  (Meyerhoff  and  Evans,  1986).  To  suppress  the  noise  due  to 
multipath  in  the  pseudorange  data,  they  are  smoothed  by  the  delta  range  (Doppler) 
data.  The  next  set  of  Equations  (3.43,  3.44,  and  3.45)  describes  this  smoothing: 

BN(tt)  =  R(tt)  -  DR(tt, ttfst)  (3.43) 


BaV8(u)  =  iBjttO/N  '  (3.44) 

Rsm(tt)  =  Bavg(“)  +  DR(tt,ufst)  (3.45) 

where, 

tt  =  GPS  time-tag  for  pseudorangc  and  delta  range  (Doppler)  data 

ttfst  =  GPS  time-tag  for  the  first  good  pseudorange  and 
Doppler  count  of  this  data  span 

R(tt)  =  Ionospheric  corrected  pseudorange  in  [km] 

DR(tt, ttfst)  =  Ionospheric  corrected  delta  range  (Doppler)  data 

N  =  number  of  data  points  in  the  span 

estimate  of  pseudorange  at  time  "ttfst" 
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BaVg  =  average  value  of  the  first  pseudorange 

RSm(u)  =  smooth  pseudorange,  if  the  absolute  value  of 

the  difference  between  the  raw  pseudorange  and 
smoothed  pseudorange  exceeds  20  [m],  the 
data  point  is  considered  bad. 

Both  the  pseudorange  and  delta  range  (Doppler)  data  must  be  corrected  for  Time 
Epoch,  clock  errors,  earth  rotation,  relativistic  effects,  and  Tropospheric  refraction. 
These  formulae  are  the  same  ones  used  for  the  iterative  least  squares  solution  in  the 
SOLTOivI  program  (Equations  3.12,  3.13,  3.14,  3.31,  3.32,  3.33,  3.34,  3.35,  3.29,  3.30). 
The  corrected  pseudorange  is  given  by, 

RC(tt)i  =  Rsm(tt)i  +  Cd(tt)j  -  ReKtOj  -  (TR(tt)j  *  (1  +  CR)]  +  Cer  (3.46) 

RDC(tt)  is  calculated  in  a  similar  manner  as  (3.46)  for  the  delta  range  (Doppler)  data. 

The  following  Equations  (3.47  to  3.51)  describe  computed  ranges  and  residuals  for 
pseudorange  and  delta  range  (Doppler)  data  (Meyerhoff,  1985):  . 

Computed  Ranges 

CR(«)i  =  ((X(t)j  -  Xr)2  +  (Y(t)j-Yr)2  +  (Z(t);  -  Zr)2]1/2  (3.47) 

for  the  range  at  the  beginning  of  the  Doppler  interval, 

CR(tt  -  td)j  =  [{X(t2)j  -  [Xr  -  (XVr  *  td)]}2  +  {Y^  -  [Yr  -  (YVr  *  td)]}2  (3.48) 

+  (Z(t2);  -  [Zr  -  (ZVr  *  td)]}  2){l2 

position  of  i^  satellite  in  an  earth- 
fixed  geocentric  cartesian  coordinate 
system  at  transit  time  t 

Xr,  Yr,  Zr  =  estimate  of  the  receiver's  position  in  an 

earth-fixed  geocentric  cartesian  coordinate 
system 


where, 

X(t),  Y(t),  Z(t)  = 
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X(t2),  Y(t2),  Z(t2)  =  position  ofi1*1  satellite  in  an 

earth-fixed  geocentric 
cartesian  coordinate  system  at 
time  t2 

XVr,  YVr,  ZVr  =  estimate  of  receiver's  velocity 

tt  =  GPS  time-tag  corrected  for  receiver  clock  bias 

td  '=  Doppler  count  interval  for  delta  range  (Doppler) 
data 

Therefore  the  computed  range  difference  for  the  individual  satellites 
CRD(tt)j  =  CR(tt)j  -  CR(tt  -  td)j 


Residuals 

The  pseudorange  residuals  are  given  by, 

Romc(tt)j  =  RC(tt)j  -  CR(tt)j  -  RCB 
and  for  the  delta  range  (Doppler)  data, 

RDomc(tt)i  =  RDC(tt)i  -  CRD(tt)j  -  RCD 
where, 

RC(tt)j  =  corrected  pseudorange  (Equation  3.46) 

RDC(tt)j  =  corrected  delta  range  (Doppler) 

CR(tt)j  =  computed  pscudorange 

DCR(tt)j  =  computed  delta  range  (Doppler) 

RCB  =  estimate  of  the  receiver  dock  bias 

RCD  =  estimate  of  frequency  bias  of  the  receiver  oscillator 
which  is  also  receiver  clock  drift. 


tt  =  GPS  time-tag  corrected  for  receiver  clock  bias 

td  =  Doppler  count  interval  for  delta  range  (Doppler) 
data 

i  =  i1*1  satellite  being  tracked 

A  state  vector  SP  is  corrected  and  improved  over  time,  is  composed  of 
the  receiver  position  in  an  earth-fixed  geocentric  cartesian  coordinate  system  (Xr,  Yr, 
Zr)  and  receiver  clock  bias  (RCB).  The  last  set  of  Equations  (3.52  to  3.58)  describes 
the  observation  equation  and  development  of  the  state  vector  for  the  pseudorange  data 
(Meyerholf,  1985): 

Observation  Equation 

D(tt)j  =  CR(tt)j  +  RCB  (3.52) 

or  rewriting  it  as, 

D(tt)i  =  [(X(t):  -  Xr)2  +  (Y(t)j  -  Yr)2  +  (Z(t)j  -  Zr)2]1^  -f-  RCB  (3.53  * 

The  H-matrix  is  defined  as  containing  the  partial  derivatives  of  the 
observation  equation  with  respect  to  the  state  vector: 

Hj  i  =  dD(tt)/3Xr  =  -(X(t)j  -  Xr)/CR(tt)j  (3.54) 

H2’i  =  0D(tt)/0Yr  -  -(Y(t)j  -  Yr)/CR(tt)j 
I-Ij’j  =  5D(tt)/dZr  =  -(Z(t)j  -  Zr)/CR(tt)j 
H4’i  =  3D(tt)/5RCB  =  1 

The  PA-matrix  is  defined  as  a  4x4  variance-covariance  matrix  for 
position  and  receiver  clock  bias,  with  off  diagonal  elements  zeroed  out  and  the  diagonal 
elements  set  to  the  initial  variance  of  Xr,  Yr,  Zr,  and  RCB,  The  Kalman  vector  K  is 
given  by  the  matrix  equation, 

K(tt)j  =  PA(tt-)  *  HTj  *  l  Hi  *  PA(lt-)  *  llTi  +  Ri  F1  (3.55) 


where, 
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PA(tt-)  =  variance-covariance  matrix  at  time  tt  before  it  is . 
updated  for  this  data  point 

R  =  variance  on  the  pseudorange  observation 

i  =  i^  satellite  being  tracked 

Once  the  PA-matrix  is  updated  it  is  defined  by, 

PA(tt+)  =  PA(tt-)  -  K(tt)j  *  Hj  *  PA(tt-)  (3.56) 


State  Vector 

The  final  two  equations  define  the  updated  state  vector  SP(tt  +  ): 

SP(tt+)  =  SP(tt-)  +  K(tt)j  *  Romc(tt)j  (3.57) 

or 

Xr(+)  =  Xr(-)  +  Kj(tt)j  *  Romc(tt)}  (3.58) 

Yr(  +  )  =  Yr(-)  +  K2(tt)i  *  Romc(tt)j 
Zr(+)  =  Zr(-)  +  K3(tt)i  *  Romc(tt)j 
RCB(  + )  =  RCB(-)  +  K4(tt)i  *  Romc(tt)j 

where  (  +  )  means  the  updated  state  vector  and  (-)  means  before  the  updated  state 
vector. 

c.  Input/Output 

There  are  two  input  files  for  the  KALMN2  program  (Figure  3.1).  The 
first  file  is  an  output  file  from  the  CON9TR  program  which  contains  records  that  arc 
aligned  with  common  time-tagged  GPS  data  items.  The  second  file  is  a  user  input  file 
which  sets  the  types  of  solutions  and/or  constraints  on  that  solution,  along  with  which 
data  solutions  will  be  printed  out.  The  user  input  guide  for  this  file  is  in  (Appendix  G). 

There  arc  three  possible  output  files  for  KALMN2  (Figure  3.1).  The 
first  file  is  the  updated  receiver  position  and  time  and  any  other  print  options  that  are 
set  by  the  user  input  file.  The  second  file  is  a  position  history  file  that  is  used  in 
printing  out  the  updated  receiver  position  and  time  file.  The  third  file  is  an  optional 
GDOP  file,  that  can  be  set  in  the  user  input  file. 
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d.  Program  Modifications  at  NPS 

The  source  code  from  NSWC  for  the  KALMN2  program  was  written  in 
FORTRAN  IV  which  has  been  modified  to  run  at  NPS  on  the  IBM-3033  mainframe 
computer.  Similar  changes  to  the  source  code  for  the  CON9TR  and  SOLTOM 
program  were  also  made  to  the  KALMN2  program  including  the  executive  files 
(Chapter  III.A.2.).  The  original  code  was  written  to  handle  both  the  precise  as  well  as 
the  onboard  broadcast  ephemeris  data  in  the  navigation  data  message  from  the  SV. 
The  test  data  for  KALMN2  was  generated  at  NSWC  using  the  broadcast  ephemeris 
only,  and  as  such,  the  modified  program  version  at  NPS  has  the  precise  ephemeris 
subroutines  and  corresponding,  call,  statements  commented  out. 

e.  Solution  Techniques  Through  Kalman  Filter 

Using  a  UDU'  factorization  of  an  eight-state  extended  Kalman  filter 
(with  smoothed  pseudorange  measurements)  in  program  KALMN2,  there  are  nine 
types  of  solutions  that  may  be  calculated  depending  on  the  user  input  file  (Appendix 
G). 

Note  that,  as  was  stated  previously  in  Chapter  III.B.2.a.,  the  velocity 
solutions  from  the  four-velocity  states  are  still  being  tested  at  NSWC. 

The  initial  receiver  position  (X,  Y,  Z)  and  clock  biases  are  updated  using 
the  corrected  state  vector  from  the  Kalman  filter  (Equations  3.57  and  3.58).  If  for  a 
given  time  line  the  absolute  value  for  the  ionospheric  correction  CIO(tt)  exceeds  20  m, 
that  data  is  considered  bad  and  flagged  as  such,  with  a  quality  vector  set  to  15.  Height 
constrained  (especially  in  marine  areas)  allows  the  user's  receiver  to  track  four  satellites 
with  poor  GDOP  and  still  be  able  to  produce  a  good  solution. 


42 


IV.  TEST  RUN  RESULTS  AND  VALIDATION 


The  testing  of  the  three  programs  CON9TR,  SOLTOM,  KALMN2  was 
accomplished  by  using  two  independent  data  sots  in  the  SEF  on  9-track  magnetic  tape 
as  received  from  NSWC.  Included  also  were  hardcopy  printouts,  listing  the  solutions 
from  those  data  sets  using  the  CON9TR,  SOLTOM,  and  KALMN2  programs  on 
NSWC's  CDC  Cyber-865  computer.  A  direct  comparison  of  the  output  listing  files 
containing  the  receiver  position  solutions  and  times  for  the  SOLTOM  and  KALMN2 
programs  from  both  the  NSWC's  CDC  Cyber-865  and  the  NPS  IBM-3033  computers 
could  then  be  made. 

After  the  original  source  codes  and  some  algorithms  were  modified  for 
compatibility  on  the  NPS  IBM-3033,  the  programs  were  test  run,  using  the  FORTVS 
compiler  with  an  auto-double  precision  option  discussed  previously  (Chapter  1 1 1. A. 2.). 

A.  STATIC  POSITION  SOLUTION  RESULTS 

Table  II  shows  the  differences  between  results  computed  at  NPS  and  NSWC 
using  the  NSWC  test  data  set  (from  Henderson  Point,  85006)  for  the  SOLTOM 
program.  Running  the  CON9TR  and  SOLTOM  programs,  position  solution  results 
(X,  Y,  Z)  were  duplicated  to  the  millimeter  level  for  solution  type  seven;  to  the 
centimeter  level  for  solution  types  one  and  two;  to  the  decimeter  level  for  solution 
types  four  and  six;  and  to  the  meter  level  for  solution  type  nine  (Appendix  E). 

The  above  differences  in  the  SOLTOM  solutions  are  considered  insignificant  due 
to  the  limitations  of  the  observed  data  (Meyerhoff,  1986). 

B.  DYNAMIC  POSITION  SOLUTION  RESULTS 

Table  III  shows  the  differences  between  results  (beginning  and  ending  solutions 
using  the  test  data  set)  computed  at  NPS  and  NSWC  using  the  NSWC  test  data  (from 
Henderson  Point,  95005)  for  the  KALMN2  program.  Running  the  CON9TR  and 
KALMN2  programs,  position  solutions  and  clock  biases  were  duplicated  to  the 
millimeter  and  10'^  second  level  respectively. 

It  should  be  noted  that  not  all  options  (Chapter  II1.B.2  and  Appendix  G)  within 
the  programs  have  been  checked  and  tested.  At  this  point  it  can  be  said  that  these 
programs  have  been  established  and  validated  on  the  IBM-3033  at  NPS  to  the  degree 
of  being  able  duplicate  NSWC's  results  if  the  same  user  options  are  followed. 
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TABLE  II 

C0N9TR  AND  SOLTOM  SOLUTION  RESULTS 


Difference  Between  NPS  and  NSWC  Solutions 


X,Y,Z 

Solution  ,  Position  N  C 
Type  (kilometers)  ( 

J ock  Bias 
.seconds) 

Clock  Drift 
(seconds) 

Number 

Of  Data  Points 

ONE  _c 

Range  Data  10  3 

Receiver  Clock 

10“9 

10"11 

1133 

TWO  _r 
Doppler  Data  10  3 
Receiver  Clock 

N/A 

10-n 

1125 

FOUR  , 

Common  Range  10  ^ 
Receiver  Clock 

N/A 

N/A 

832 

SIX 

Range  Data  10  H 

S V  Clock 

10"7 

10-12 

1133 

SEVEN  _7 

Doppler  Data  10  ' 

SV  Clock 

N/A 

10-n 

1125 

NINE  o 

Common  Range  10 

SV  Clock 

10"6 

IQ-12 

832 
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TABLE  III 

C0N9TR  AND  KALMN2  SOLUTION  RESULTS 


Difference  Between  NPS  and  NSWC  Solutions 


GPS  Time-Tag 
(seconds) 

(lilome^ers) 

Clock  Bias 
(seconds) 

Number  of 
Doppler  Edited 
Points 

First  .Solution 

14129.00 

10"5 

10"9 

N/A 

Last  Solution 

17499. 00 

10“5 

h-» 

O 

1 

CT> 

9 

V.  SUMMARY 


How  to  process  the  GPS  satellite  data  from  raw  cassettes  off  the  TI-4100 
GEOSTAR  receiver  to  a  reduced  position  and  time  solution  has  been  discussed. 

The  capability  of  processing  raw  data  from  cassettes  to  SEF  at  NPS  is  a  task  for 
the  future. 

NPS  at  this  time  has  the  capability,  starting  with  a  SEF  9-track  magnetic  tape,  to 
use  the  CON9TR  program  (where  the  common  time-tagged  data  are  aligned),  to  select 
either  the  SOLTOM  (an  iterative  least  squares  for  positioning  a  static  receiver)  or  the 
KALMN2  program  (which  uses  a  Kalman  filter  for  positioning  a  dynamic  receiver)  and 
to  compute  position  and  time  solutions.  The  calculations  for  program  CON9TR,  and 
the  math  models  for  programs  SOLTOM  and  KALMN2  were  documented.  The 
changes  required  to  modify  the  original  codes  and  algorithms  for  these  programs 
written  for  NSWC's  CDC  Cyber-865  computer  to  make  them  compatible  with  the 
IBM-3033  computer  at  NPS  have  been  described.  Testing  of  the  programs  shows  a 
successful  conversion  (with  the  restrictions  discussed  in  Chapter  IV.)  of  the  programs 
to  the  IBM-3033  mainframe  computer.  This  software  forms  a  basis  for  a  GPS 
programs  library  that  will  aid  research  using  the  GPS  at  NPS. 
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VI.  RECOMMENDATIONS 


Since  the  conversion  of  the  CON9TR,  SOLTOM,  and  KALMN2  programs  has 
been  successfully  tested  using  the  two  independent  data  sets  from  NSWC,  it  is 
recommended  to  explore  other  options  within  the  programs  set  by  the  user  input  files 
for  the  SOLTOM  and  KALMN2  programs.  Options  that  should  be  explored  for  the 
SOLTOM  program  include  the  use  of  the  precise  ephcmcris  and  clock  files,  the  use  of 
more  than  one  iteration  for  the  solutions,  and  the  calculation  of  Newton- Raphson 
range  solutions  followed  by  a  comparison  with  the  iterative  batch  least  squares 
solutions  for  a  given  time  line.  The  KALMN2  program  should  be  explored  by  using 
the  precise  ephemeris  and  clock  files,  turning  the  height  constraint  on  and  off  in  the 
user  input  file  and  comparing  the  solution  results.  The  refraction  correction  should  be 
turned  on  and  off  to  sec  how  the  solutions  are  affected  (this  is  another  option  in  the 
user  input  file). 

It  is  further  recommended  to  explore  the  possiblities  of  rewriting  or  adding  codes 
and  algorithms  to  improve  the  solution  corrections  for  updating  the  four-state  vectors 
of  the  velocity  in  the  extended  eight-state  Kalman  filter  making  them  reliable  for  use  in 
velocity  solutions. 

In  the  future  TI-4100  software  for  carrier  phase  data,  differential  positioning  etc., 
should  be  acquired  and  modified  to  be  compatible  with  the  IBM-3033  computer  at 
NPS. 
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APPENDIX  A 

C0N9TR  SUBROUTINE  DESCRIPTIONS  AND  NESTING 


INPUT  - 
GETDAT - 

RDDTL  - 

ER9TR  - 
UDNTDT - 

PRTSL  - 
UDWTDT  - 

UDSATD  - 
FDTRN  - 
DEDPOF- 

PRTPDT  - 
UDTRD  - 

DEQUDT - 


DERGDT - 

DEDPDT - 


This  routine  reads  the  user  inputs. 

This  routine  reads  the  next  set  of  satellite 
tracking  data  from  the  9-track  Standard  tape 
and  also  updates  weather,  Navigation  data, 
almanac  and  calibration. 

This  routine  initializes  the  quality  vector 
llags  and  controls  branching  of  the  program 
depending  on  a  particular  error  Hag. 

This  routine  controls  branching  of  the  program 
for  several  error  and  terminating  conditions. 

This  routine  processes  non-tracking  data  so  that 
it  can  be  used,  after  it  determines  which  type 
of  non-tracking  data  was  read  from  the  9-track 
Standard  tape. 

This  routine  prints  tape  solutions. 

This  routine  acts  the  weather  data  from  the  9-track 
Standard  tape  and  stores  it  in  common  that  will  be 
used  to  calculate  the  tropospheric  correction. 

This  routine  'eads  the  9-track  Standard  tape  and 
defines  new  satellite  orbit  clock  parameters. 

This  routine  will  give  the  tracker  number  for  the 
satellite  that  the  new  data  is  for. 

This  routine  stores  the  values  for  the  LI  and  L2 
Doppler  offsets  so  that  they  may  be  used  in  the 
program. 

This  routine  prints  any  unused  9-track  Standard 
tape  data. 

This  routine  decides  if  the  tracking  data  is 
continuous  count  Doppler  range  or  data  quality 
measurement  blocks  and  then  calls  the  correct 
subroutine  to  handle  the  correct  data  type. 

This  routine  defines  data  qualitv  and  then  stores 
the  data  quality  information.  The  satellite 
tracker  mode  and  the  quality  vectors  arc  used  to 
determine  if  the  data  are  good.  The  remaining  data 
values  in  this  block  arc  used  to  calculate  errors 
on  the  tracking  data. 

This  routine  stores  the  range  data  and  the  signal 
to  noise  data  in  correct  data  array  that  will  Ue 
processed. 

This  routine  converts  the  accumulated  Doppler 

Rhase  data  into  range  difference  data.  The 
lopplcr  phase  data  is  corrected  for  offset  and 
then  it  is  dilfcrcnccd  with  preceding  data  if  the 
preceding  data  was  good. 
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GNINP  - 


This  routine  reads  the  logical  unit  number  of  the 
input  file  and  returns  the  information  read  one 
item  at  a  time. 

GNID  -  This  routine  controls  how  data  items  are  read  by 

reading  a  four  digit  data  type  number. 

IN9TR  -  This  routine  initializes  the  read  package  of  the 

9-track  Standard  tape  and  sets  up  a  data  block 
for  the  data  items  read  off  of  the  tape. 

BLDATA  -  This  routine  contains  the  data  that  has  been  read 
off  of  the  9-track  Standard  tape. 

SDAT  -  This  routine  takes  the  site  name  and  the  data 

interval  from  the  9-track  Standard  tape  read 
common  and  passes  these  values. 

FSIG  -  This  routine  determines  the  difference  between 

the  raw  ranae  and  the  range  difference  data  for 
the  satellite  oeing  tracked. 

BLOCK  DATA  -  This  routine  was  added  at  NPS  to  make  this 
program  compatible  with  the  FORTRAN  77 
compilers  on  the  IBM-3033  mainframe  computer. 
It  consists  of  common  statements  and  data 
statements  that  the  rigorous  compilers 
at  NPS  would  not  permit  where  they  were  placed 
in  the  original  source  code  of  CON9TR. 


MAIN  PROGRAM  CON9TR 
INPUT 
GETDAT 

RDDTL 

ER9TR 

UDNT.DT 

PRTSL 

UDWTDT 

UDSATD 

FDTRN 

DEDPOF 

PRTPDT 

UDTRD 

DEQUDT 

FDTRN 

DERGDT 

FDTRN 

DEDPDT 

FDTRN 

GNINP 

GNID 

IN9TR 

BLDATA 

GNINP 

RDDTL 

SDAT 


FSIG 

BLOCK  DATA 
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APPENDIX  B 

INPUT  GUIDE  FOR  PROGRAM  CON9TR 


(  by  S.  L.  MeycrhofT,  NSWC) 


This  is  the  user  input  file  for  the  program  with  variable  names  and  descriptions.  The 
input  is  format  free. 


line  I 

1:  NSAT 


number  of  trackers  (4) 


2:  TSS 


GPS  time  of  first  data  point  to  put  on 
new  file 


3:  TES 


4:  K2 
5:  K4 
line  2 


1:  X 
2:  Y 
3:  Z 


time  to  end  data  on  new  file,  if  EOF 
found  before  TES  then  data  will  end  at 
EOF 

range  variance  bias  factor  (kilometer}^ 
doppler  variance  bias  factor  (kilometer)^ 


X  coordinate  of  initial  receiver  position 
that  will  be  put  on  the  new  files 

Y  coordinate  of  initial  receiver  position 
that  will  be  put  on  the  new  files 

Z  coordinate  of  initial  receiver  position 
that  will  be  put  on  the  new  files 
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APPENDIX  C 

TEXAS  INSTRUMENTS  GPS  SATELLITE  RECEIVER  FILE  FORMAT 

(  by  S.  L.  Meyerhoff,  1984,  NSWC  ) 

Data  from  the  DMA  version  of  the  TI-4100  GPS  satellite  receiver  can  be  put 
into  the  following  format.  In  this  format  all  common  time  satellite  tracking  data  will 
be  in  the  same  record  block.  The  format  for  this  file  is  as  follows. 


var 

type 

meaning 

Record  One 

word  1  ITYPE 

intg 

This  is  the  data  type  for 

the  next  record 

Record  Two 

If  ITYPE  equals  one  then  record  two  will  contain  initial  information. 


word  1  RL10FF 

real 

frequency  bias  ofiset  for  L 1' 

word  2  RL20FF 

real 

frequency  bias  ofiset  for  L2 

word  3  PB 

real 

pressure  in  (millibars) 

word  4  TP 

real 

temperature  (degrees  Celsius) 

word  5  HD 

real 

relative  humidity  (percent) 

word  6  X 

real 

X  coordinate  of  initial 

receiver  position  (kilometers) 

word  7  Y 

real 

Y  coordinate  of  initial 

receiver  position  (kilometers) 

word  8  Z 

real 

Z  coordinate  of  initial 

receiver  position  (kilometers) 

word  9  NTR 

intg 

maximum  number  of  trackers  used 

will  be  4  or  larger 

word  10  DTREC 

real 

user  given  data  interval 

(seconds) 

word  11  SITE(l) 

char 

site  identification  for  data 

word  12  SITE(2) 

char 

site  identification  for  data 
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(continued) 


If  ITYPE  equals  two  then  record  two  will  contain  satellite  ephemeris  data  (navigation 
data  message) 


word  1  J 

intg 

tracker  number 

word  2  ISAT(J) 

intg 

satellite  prn  number 

word  3  AS(1) 

real 

satellite  clock  bias  in  seconds 

word  4  AS(2) 

real 

satellite  clock  drift 

(seconds/second) 

word  5  AS(3) 

real 

satellite  clock  aging 
(seconds/second*second) 

word  6  CICS 

real 

cross  track  cosine  amplitude 
(radians) 

word  7  CISS 

real 

cross  track  sine  amplitude 
(radians) 

word  8  CRCS 

real 

radial  cosine  amplitude 
(kilometers) 

word  9  CRSS 

real 

radial  sine  amplitude 
(kilometers) 

word  10  CUCS 

real 

in  track  cosine  amplitude 
(radians) 

word  11  CUSS 

real 

in  track  sine  amplitude 
(radians) 

word  12  DNS 

real 

corrections  to  mean  motion 

(radians) 

word  13  ES 

real 

eccentricity 

word  14  I DOTS 

real 

inclination  rate 

(radians/second) 

word  15  IOS 

real 

inclination  (radians) 

word  16  MOS 

real 

mean  anomaly  (radians) 

word  17  OMEDS 

real 

right  ascension  rate 
(radians/sccond) 

word  18  OMEGS 

real 

right  ascension  (radians) 

word  19  SQAS 

real 

square  root  of  semi-major 
axis  (kilometers) 
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word  20  ADE 

real 

age  of  data  ephemeris 
(seconds) 

word  21  TOCS 

real 

clock  epoch  (seconds) 

word  22  TOES 

real 

ephemeris  epoch  (seconds) 

word  23  IWKNOS 

intg 

week  number 

word  24  WS 

real 

argument  of  perigee 
(radians) 

word  25  IEDATS 

intg 

sv  health 

word  26  ADC 

real 

age  of  data  clock  (seconds) 

If  ITYPE  equals  three  then  record  two  will  contain  the  tracking  information. 


word  1  TT 

real 

GPS  time-tag  for  data 

that  follows  (seconds) 

word  2  ISAT(l) 

The  PRN  number  for  the 

to 

to 

intg 

satellite  being  tracked 

word  NTR+  1  ISAT(NTR) 

on  trackers  i  to  NTR 

word  NTR+2 

CR1(1) 

The  LI  pseudorange  on 

to 

to 

real 

each  tracker  (seconds) 

.word  NTR*2+  1 

CRl(NTR) 

word  NTR*2  +  2  CR2(1) 

The  L2  pseudorange  on 

to 

to 

real 

each  tracker 

word  NTR*3  + 1 

CR2(NTR) 

(seconds) 

word  NTR*3  +  2 

DOPl(l) 

The  LI  continuous 

to 

to 

real 

count  dopplcr  data  on 

word  NTR*4  +  1 

DOPl(NTR) 

each  tracker  (cycles) 

word  NTR*4  +  2 

DOP2(l) 

The  L2  continuous 

to 

to 

real 

count  dopplcr  data  on 

word  NTR*5+  1 

DOP2(NTR) 

each  tracker  (cycles) 

word  NTR*5  +  2 

SGR1(1) 

The  variance  on  the  LI 

to 

to 

real 

pscudorange 

word  NTR*6+  1 

SGRl(NTR) 

(kilometcr’kilometcr) 

word  NTR*6+2 

SGR2(1) 

The  variance  on  the  L2 
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to  to 

word  NTR*7  +  1  SGR2(NTR) 

word  NTR*7 + 2  SGDl(NTR) 
to  to 

word  NTR*8+  1  SGDl(NTR) 

word  NTR*8  +  2  SGD2(1) 
to  to 

word  NTR*9  +  1  SGD2(NTR) 

word  NTR*9  +  2  MQVEL(l) 
to  to 

word  NTR*10+  1  MQVEL(NTR) 


real  pseudorange 

(kilometer*kilometer) 

The  variance  on  the  LI 

real  doppler 

(kilometer*kilometer) 

The  variance  on  the  L2 

real  doppler 

(kilometer*kilometer) 

The  quality  vector  for 
intg  all  trackers 

15  =  loss  of  lock 
14  =  bad  data,  no 
loss  of  lock 
1  =  good  range  data 
bad  doppler  data 
0  =  all  data  good 
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APPENDIX  D 

SOLTOM  SUBROUTINE  DESCRIPTIONS  AND  NESTING 


CARDIO  - 
LLXYZ - 

SIGXYZ  - 

MATRIX  - 
LINV2F  - 
VMULFF  - 
VMULFM  - 
INDAT  - 
NEWSOL  - 
LSTSOL  - 

GETDAT - 

DEDPDT - 

DEFBR  - 
CRD AT  - 

BSADR  - 

CRDAT2  - 
ION  - 
IONRD  - 


This  routine  reads  the  input  for  program  SOLTOM 
from  the  output  of  CONyTR. 

This  routine  converts  position  in  the  geodetic 
coordinate  system  to  tne  earth-fixed  geocentric 
coordinate  system. 

This  routine  converts  sigmas  in  the  geodetic 
system  into  sigmas  in  the  earth-fixed 
geocentric  coordinate  system. 

This  routine  performs  matrix  operations  by 
calling  routines  in  the  IMSL  library. 

This  routine  is  in  the  IMSL  and  does  a  matrix 
inversion. 

This  routine  is  in  the  IMSL  and  does  matrix 
multiplication. 

This  routine  is  in  the  IMSL  and  does  a 
transpose  of  matrix  multiplication. 

This  routine  will  pass  initial  values  to  some 
data  arrays. 

This  routine  sets  up  the  configuration  for  a 
for  a  given  solution  routine. 

This  routine  calculates  a  sequential  least 
squares  solution  for  the  receiver  position  and 
clock. 

This  routine  reads  the  next  set  of  satellite 
tracking  data  from  CON9  i'R  output  file  and  it 
will  also  update  weather,  Navigation  data, 
almanac  and  calibration. 

This  routine  converts  the  accumulated  Doppler 
phase  data  into  range  difference  data  by 
correcting  phase  for  offset  and  then 
differencing  it  with  data  from  before 
(if  before  data  is  good). 

This  routine  defines  bias  range 

This  routine  defines  correction  for  receiver 
bias. 

This  routine  will  correct  data  for  receiver 
clock  bias,  drift  and  also  correct  GPS  time 
tag  for  the  end  of  week  cross  over. 

This  routine  converts  data  into  range 
residuals  and  corrects  data  for  error  sources. 

This  routine  corrects  for  ionospheric 
correction. 

This  routine  corrects  range  difference  data. 
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EPI-IDAT  - 


CLOCK  - 

SATPOS - 

ESTPOS  - 
ZANG  - 

RES  - 

SMOT  - 
TROP  - 
REL  - 
FSIG  - 

SOLTNS - 

SMPT1  - 

UDRRC  - 

SMPT2  - 

UDRDRC - 

SMPT3  - 


This  routine  defines  the  following  values 
from  the  ephemeris  data,  satellite  tracker 
number,  satellite  bias,  drift  and  aging  terms, 
semi-major  axis  of  satellite  orbit, 
eccentricity  of  the  satellite  orbit  and  epoch 
time  for  clock  terms. 

This  routine  adds  satellite  and  clock  error 
and  converts  it  to  range  measurements  in 
kilometers  and  also  performs  the  same  for  range 
difference  data. 

This  routine  will  find  the  position  of  the 
satellite  from  the  Navigation  Message  or  from 
the  Precise  Ephemeris. 

This  routine  finds  the  position  of  the 
satellite  from  the  Precise  Ephemeris  data. 

This  routine  finds  the  Zenith  angle  by  doting 
the  X,  Y,  Z  positions  of  the  satellite  receiver 
position  vector  to  the  vector  from  the  receiver 
to  the  satellite. 

This  routine  finds  the  residual  by  calculating 
the  range  to  a  satellite  and  subtracting  it 
from  the  observed  ranae.  This  is  also 
performed  for  range  difference  data. 

This  routine  determines  an  earth  rotation 
correction. 

This  routine  finds  tropospheric  refraction 
correction. 

This  routine  finds  a  relativistic  correction 
to  the  data. 

This  routine  finds  the  variancefsigma  squared) 
for  the  raw  range  and  range  difference  data  for 
the  satellite  being  tracked. 

This  routine  calls  the  different  routines  to 
sum  data  in  the  correct  part  of  the  "B"  matrix 
and  "E"  vector  depending  on  which  options  have 
been  set  by  the  user  for  type  of  solution. 

This  routine  calls  a  routine  to  define  additions 
to  "B"  matrix  and  "E"  vector  plus  directing  the 
called  routine  on  where  to  store  this 
information  for  the  "B"  matrix  and  "E"  vector 
(for  range  data  only). 

Update  range  solution  matrix  receiver  clocks. 
Then  store  data  in  matrix  so  that  solution  can 
computed. 

This  routine  calls  a  routine  to  define 
additions  to  "B"  matrix  and  "E"  vector  plus 
directing  the  called  routine  on  where  to  store 
this  information  for  "B"  matrix  and  "E"  vector 
(for  range  difference  data  only) 

Update  range  difference  solution  matrix 
receiver  bias.  This  routine  stores  data  in 
matrix  so  a  solution  can  be  computed. 

This  routine  calls  a  routine  to  define 
additions  to  "B"  matrix  and  "E"  vector,  plus 
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NETRAP - 

CROUT  - 
LINV3F  - 
POSAVG  - 
SMPT6  - 

UDRSC  - 

BELIM  - 
BMSTR  - 
SMPT7  - 

UDRDSC  - 

SMPT8  - 

SMPT9  - 

UDCRSC  - 

ELIM  - 
SMPTIO  - 


directing  the  called  routine  on  where  to  store 
this  inlormation  lor  the  "B"  matrix  and  "E" 
vector  (for  range  and  range  difference  data 
only). 

This  routine  uses  a  Newton-Raphson  method  to 
calculate  time  bias  adjustment  and  then  improve 
the  receiver  initial  position. 

This  routine  calls  an  IMSL  routine  that  finds 
the  determinate  of  a  matrix. 

This  is  an  IMSL  routine  that  calculates  the 
determinate  of  a  matrix. 

This  routine  finds  the  average  of  all  positions 
up  to  current  time  line. 

This  routine  calls  a  routine  to  define  additions 
to  "B"  matrix  and  "E"  vector  plus  directing  the 
called  routine  on  where  to  store  this 
information  (for  range  data  only,  satellite 
biases) 

This  routine  stores  data  in  matrix  so  solution 
so  solution  can  be  computed  after  the  range 
solution  matrix  and  satellite  clocks  have  been 
updated. 

This  routine  eliminates  certain  parameters  in 
the  "B"  matrix. 

This  routine  stores  data  in  a  matrix  so  a 
solution  can  be  computed. 

This  routine  calls  a  routine  to  define 
additions  to  the  "B"  matrix  and  "E"  vector 
plus  directing  the  called  routine  on  where  to 
store  this  information  (for  range  difference 
data  only,  satellite  clock). 

This  routine  stores  data  in  a  matrix  so  that 
a  solution  can  be  computed  alter  updating  the 
range  solution  matrix  and  satellite  bias. 

This  routine  calls  a  routine  to  define  additions 
to  "B"  matrix  and  "B"  vector  plus  directing  the 
called  routine  where  to  store  this  information 
(for  range  and  range  difference  data  only). 

This  routine  calls  a  routine  to  define  additions 
to  the  "B"  matrix  and  "E"  vector  plus  directing 
the  called  routine  on  where  to  store  this 
information  (for  common  time  range  data  only). 

This  routine  stores  data  in  a  matrix  so  that 
a  solution  can  be  computed  after  updating  the 
common  range  solution  matrix  and  satellite 
clocks. 

This  routine  finds  data  needed  for  elimination 
and  then  calls  a  routine  to  do  that  elimination. 

This  routine  calls  a  routine  to  define 
additions  to  the  "B"  matrix  and  "E"  vector  plus 
directing  the  called  routine  on  where  to  store 
this  information  (for  Doppler  data  used  as  bias 
range  data). 
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BSRNG  - 


NEW ITR  - 
BLOCK  DATA  - 


This  routine  defines  additions  to  the  "B" 
matrix  and  "E"  vector. 

This  routine  resets  data  so  that  a  new  iteration 
can  be  done. 

This  routine  was  added  at  NPS  to  make  this 
program  compatible  with  the  FORTRAN  77 
compilers  on  the  IBM-3033  mainframe 
computer  system.  It  contains  common 
statements  and  data  statements  that 
NPS's  compilers  would  not  permit  where 
they  were  placed  in  the  original  code  of 
program  SOLTOM. 


MAIN  PROGRAM  SOLTOM 
CARDIO 
LLHXYZ 
SIGXYZ 
MATRIX 

LINV2F 

VMULFF 

VMULFM 

INDAT 

NEWSOL 

LSTSOL 

MATRIX 

GETDAT 

DEDPDT 

DEFBR 

CRD  AT 

BSADR 

CRDAT2 

ION 

IONRD  . 

EPIIDAT 

CLOCK 

SATPOS 

ESTPOS 

ZANG 

RES 

SMOT 

'PROP 

REL 

FSIG 

SOLTNS 

SMPTI 

UDRRC 

SMPT2 

UDRDRC 

SMPT3 

UDRRC 

UDRDRC 

NETRAP 

CROUT 

LINV3F 

CRDAT 

POSAVG 

SMPT6 

UDRSC 

BELIM 

MATRIX 

LINV2F 

VMULFF 

VMULFM 

BMSTR 

SMPT7 


5S 


UDRDSC 

BELIM 

BMSTR 

SMPT8 

UDRSC 

UDRDSC 

SMPT9 

UDCRSC 

ELIM 

BELIM 

BMSTR 

SMPT10 

BSRNG 

CRDAT 

UDRSC 

NEWSOL 

NEWITR 

GETDAT 
BLOCK  DATA 
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APPENDIX  E 

INPUT  GUIDE  AND  SOLUTION  TYPES  FOR  PROGRAM  SOLTOM 


(  by  S.  L.  Mcycrhoff,  NSWC) 

INPUT  GUIDE 

Var  Type  Description 

CARD  ONE  Format  (2F 10. 0,615) 

(This  line  can  be  read  format  free) 

col  1  to  10 

TSS  real  The  GPS  time  tag  for  data  at  start 

of  solution  (seconds) 

col  1 1  to  20 

TES  real  The  GPS  time-tag  of  last  data  point 

to  use  in  so'ution  (seconds) 

col  21  to  25  (right  justified) 

NPN  integer  The  number  of  time  lines  of  data 

needed  before  first  solution 

col  26  to  30  (right  justified) 

NPli  integer  The  number  of  time  lines  of  data 

between  solutions 

col  31  to  35  (right  justified) 

ITRY  integer  The  maximum  number  of  iterations 

col  36  to  40  (right  justified) 

NSAT  integer  The  maximum  number  of  trackers  to 

process  at  once  (range  1  to  4) 

col  41  to  45  (right  justified) 

IG  f  integer  =  1  if  data  source  is  Standard  9-track 

format 

=  2  if  data  source  GESAR  test  file 
format 

=  3  if  data  source  is  NAV  IF  tape 

col  46  to  50  (right  justified) 

ISP  integer  =  1  for  satellite  position  from  NAV 

data 

=  2  for  satellite  position  from 
Celestial  Earth  fixed  trajectory 
files 

col  51  to  55  (right  justified) 

IAGE  integer  =  1  solve  for  clock  aging,  otherwise 

not  solved  for 

col  56  to  60  (right  justified) 

ICLOCK  integer  =  0  solve  clock  corrections  are  from 

NAV  data  message 
=  1  clock  corrections  are  post  fit 
data  from  cards 

CARD  TV/0  ( receiver  position )  Formal  ( 1 10  JG  10.0 ) 

('fins  card  can  be  read  format  free) 

col  10 
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I  FLAG 

integer 

=  0  for  cartesian  coordinates 
=  1  for  geodetic  coordinates 

col  11  to  20 

XL  AT 

real 

latitude  or  X  coordinate  for  position 
(degrees  or  kilometers) 

col  21  to  30 
YLOG 

real 

longitude  or  Y  coordinate  for  position 
(degrees  or  kilometers) 

col  31  to  40 
ZHT 

real 

height  or  Z  coordinate  for  position 
(kilometers) 

col  41  to  50 

TB 

real 

time  bias  for  receiver  clock 
(kilometers) 

col  51  to  60 

TD 

real 

time  drift  for  receiver  clock 
(kilometers/second) 

col  61  to  70 

TO 

real 

epoch  for  clock  drift  (given  and 
solution)  If  negative,  then  first 
used  time  line  is  the  epoch  for  the 
solution  clock  drift 

col  71  to  80 
AGE 

real 

aging  for  receiver  clock 
(kilometers/secondvsecond) 

CARD  THREE  ( Reference  Geoid)  Format  ( 2GI0.0 ) 

(This  card  can  be  read  format  free) 


col  1  to  10 

AG  real  semi-major  axis  of  the  reference  geoid 

(kilometers) 


col  1 1  to  20 

ROBL  real  Obliquity  of  the  reference  geoid 

CARD  FOUR  {Sigmas  for  Position )  Format  [7 G  10.0) 

(This  card  can  be  read  format  free) 


col  1  to  10 
S1GLT 

real 

sigma  of  X  or  latitude  coordinate 
(kilometers) 

col  11  to  20 
SIGLG 

real 

sigma  of  Y  or  longitude  coordinate 
(kilometers) 

col  21  to  30 
SIGHT 

real 

sigma  of  Z  or  height  coordinate 
(kilometers) 

col  31  to  40 
SIG(I) 

real 

sigma  of  refraction  scaling  constant 

col  41  to  50 
S1G(2) 

real 

sigma  on  receiver  or  satellite  1  time 
bias  (kilometers) 

col  51  to  60 
SIG(3) 

real 

sigma  on  receiver  or  satellite  1 
clock  drift  (kilometcrs/sccond) 
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real 


col  61  to  70 


SIG(4) 


sigma  for  aging  term  or  satellite  2 
bias  (kilometers/sccond*sccond) 


CARD  FIVE  ( Sigmas  for  Position )  Format  (7GI0.0) 

(This  card  can  be  read  format  free) 


col  1  to  10 
SIG(5) 

real 

siema  on  satellite  2  drift 
(kilometers/second) 

col  1 1  to  20 
SIG(6) 

real 

siema  for  satellite  3  bias 
(kilometers) 

col  21  to  30 
SIG(7) 

real 

siema  for  satellite  3  drift 
(Kilometers/second) 

col  31  to  40 
SIG(8) 

real 

sigma  for  satellite  4  bias 
(kilometers) 

col  41  to  50 
SIG(9) 

real 

sigma  for  satellite  4  drift 
(kilometers) 

SOLUTION  TYPES 

CARD  SIX  ( Types  of  Solutions  to  do)  Format  (1615) 

(This  card  cannot  be  read  format  free) 

col  5 

NSOL(l) 

integer 

=  1  for  range  solution  receiver  clock 

0  =  no 

col  -10 

NSOL(2) 

integer 

=  1  for  range  difference  solution 
receiver  clock 
=  0  no 

col  15 

NSOL(3) 

integer 

=  1  for  range  and  range  difference 
solution  receiver  clock 
=  0  no 

col  20 

NSOL(4) 

integer 

=  1  for  common  time  range  solution 

receiver  clock 

=  0  no 

col  25 

NSOL(5) 

integer 

=  1  for  Newton-Raphson  range  solution 
=  0  no 

col  30 

NSOL(6) 

integer 

=  1  for  range  solution  satellite 

clocks 

=  0  no 

col  35 

NSOL(7) 

integer 

=  1  for  range  difference  solution 

satellite  clock 

=  0  no 

col  40 
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NS0L(8) 

integer 

=  1  for  range  and  range  difference 

solution  satellite  clock 

=  0  no 

co!  45 

NSOL(9) 

integer 

=  1  for  common  time  range  solution, 

satellite  clock 

=  0  no 

col  50 
NSOL(IO) 

integer 

=  1  for  doppler  data  used  as  bias 

range  data 

=  0  no 

CARD  SEVEN  ( Variance  Dias  Factor)  Format  ( 2G10.0 ) 

(This  card  can  be  read  format  free) 

col  1  to  10 

K2 

real 

bias  factor  for  variance  on  range 
data 

col  11  to  20 

K4 

real 

bias  factor  for  variance  on  doppler 
data 

CARD  EIGHT  ( Print  Options )  Format  (515) 

(This  card  can  be  read  format  free) 

col  5 

IPRT(l) 

integer 

=  1  print  position  solution  (X,Y,Z) 

=  0  do  not  print 

col  10 

IPRT(2) 

integer 

=  1  print  solution  sigmas  and  clock 

terms 

=  0  do  not  print 

col  15 

IPRT(3) 

integer 

=  1  print  data  and  time  tags 
=  0  do  not  print 

col  20 

IPRT(4) 

integer 

=  1  print  residuals,  range,  trop, 

data  sigma 

=  0  do  not  print 

col  25 

IPRT(5) 

integer 

=  1  print  satellite  position 
=  0  do  not  print 

col  30 

IPRT(6) 

integer 

=  1  print  residuals  and  all  data 

corrections 

=  0  do  not  print 

col  35 

1PRT(7) 

integer 

=  1  print  solution  position  in 

lat,  long,  height 

=  0  do  not  print 

col  40 

1PRT(8) 

integer 

=  1  print  real  time  solutions  from 
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9-track 


=  0  do  not  print 


CARD  NINE  ( Weather  Data,  Elev.  Cutoff)  Format  ( 4G10.0 ) 

(This  card  can  be  read  format  free) 


col  1  to  10 

PB 

real 

pressure  in  millibars 

col  11  to  20 

TC 

real 

temperature  in  degrees  Centigrade 

col  21  to  30 

HD 

real 

relative  humidity  in  percent 

col  31  to  40 

EL 

real 

elevation  angle  cutoff  (degrees) 

CARD  TEN  ( Station  Number)  Format  (15) 

(This  card  can  be  read  format  free) 

col  1  to  5 

ISTA 

integer 

station  identification 

CARD  ELEVEN  (Tracker  1  Clock  Parameters)  Format  (3G10.0) 

col  1  to  10 

AO 

real 

clock  bias  of  satellite  tracked  on 
tracker  1 

col  1 1  to  20 

A1 

real 

clock  drift  of  satellite  tracked  on 
tracker  1 

col  21  to  30 
,  A2 

real 

clock  aging  of  satellite  tracked  on 
tracker! 

CARD  TWELVE  (Tracker  2  Clock  Parameters)  Format  (3G10.0) 

('1  his  card  can  be  read  format  free) 

col  1  to  10 

AO 

real 

clock  bias  of  satellite  tracked  on 
tracker  2 

col  1 1  to  20 

A1 

real 

clock  drift  of  satellite  tracked  on 
tracker  2 

col  21  to  30 

A2 

real 

clock  aging  of  satellite  tracked  on 
tracker! 

CARD  THIRTEEN  (Tracker  3  Clock  Parameters)  Formal  (3G10.0) 

(This  card  can  be  read  format  free) 

col  1  to  10 

AO 

real 

clock  bias  of  satellite  tracked  on 
tracker  3 

col  11  to  20 

A1 

real 

clock  drift  of  satellite  tracked  on 
tracker  3 

col  21  to  30 

A2 

real 

clock  aging  of  satellite  tracked  on 
tracker! 
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CARD  FOURTEEN  ( Tracker  4  Clock  Parameters )  Format(3G10.0 ) 

(This  card  can  be  read  format  free) 

col  1  to  10 


AO 

real 

clock  bias  of  satellite  tracked  on 
tracker  4 

col  11  to  20 

A1 

real 

clock  drift  of  satellite  tracked  on 
tracker  4 

col  21  to  30 

A2 

real 

clock  aging  of  satellite  tracked  on 
tracker  4 
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APPENDIX  F 

KALMN2  SUBROUTINE  DESCRIPTIONS  aND  NESTING 


CARDIO  - 
LLHXYZ  - 

SIGXYZ  - 

MATRIX  - 

LINV2F  - 

VMULFF  - 

VMULFM  - 

INDAT  - 

TIMEUP  - 

PROMPT  - 
GETDAT - 

NEWRAN  - 

DEDPDT - 

DEFBR  - 
CRDAT  - 

BSADR  - 

CRDAT2  - 

ION  - 

IONRD  - 
EPMDAT - 

CLOCK  - 


This  routine  reads  input  for  program. 

This  routine  converts  latitude,  longitude  and 
height  to  X,  Y,  Z  in  earth-fixed  geocentric 
cartesian  coordinates. 

This  routine  converts  sigmas  in  the  geodetic 
coordinate  system  into  sigmas  in  the  earth- 
fixed  geocentric  coordinate  system. 

This  routine  performs  matrix  operations  by 
calling  the  proper  routine  in  the  IMSL  library. 

This  routine  is  in  the  IMSL  and  performs  a 
matrix  inversion. 


This  routine  is  in  the  IMSL  and  oerforms  a 
matrix  multiplication. 


This  routine  performs  a  transpose  of 
multiplication  ol  matrices. 


This  routine  gives  the  initial  values  to  the 
"B"  matrix  (U  in  Kalman  filter) 


This  routine  uses  a  UDU'  factorization  method 
to  update  state  matrix. 


This  routine  aids  in  the  UDU'  factorization. 


This  routine  reads  the  next  set  of  satellite 
tracking  data  from  CON9TR  output  file. 

This  routine  reduces  noise  and  multipath  in  the 
range  data. 

This  routine  converts  Doppler  phase  into  range 
differenced  data. 


This  routine  defines  bias  range. 

This  routine  defines  corrections  for  receiver 
bias. 


This  routine  will  correct  data  for  receiver 
clock  bias  and  drift  and  also  correct  GPS  time- 
tag  for  the  end  of  week  cross  over. 

This  routine  converts  data  into  range  residuals 
and  corrects  data  lor  error  sources. 


This  routine  corrects  for  ionospheric 
refraction. 


This  routine  corrects  range  difference  data 

This  routine  defines  the  values  from  the 
ephemeris  data. 

This  routine  adds  satellite  and  clock  error  and 
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converts  it  to  range  measurements  in  kilometers 
and  also  perforrti  the  same  lor  range  difTcrence 
data. 

SATPOS  - 

This  routine  will  find  the  position  of  the 
satellite  from  the  navigation  data  message  or 
from  the  Celestial  Earth-Fixed  Trajectory 
(Precise  Ephcmeris). 

ESTPOS  - 

This  routine  finds  the  position  of  the 
satellite  from  the  Precise  Ephcmeris  data. 

ZANG  - 

This  routine  finds  the  zenith  angle  by  doting 
the  X,  Y,  Z  positions  of  the  satellite  receiver 
position  vector  to  the  vector  from  the  receiver 
to  the  satellite. 

RES  - 

This  routine  finds  the  residual  for  ranges. 

SMOT  - 

This  routine  determines  an  earth  rotation 
correction. 

TROP  - 

This  routine  calculates  a  tropospheric 
refraction  correction. 

REL  - 

This  routine  finds  a  relath  istic  correction 
to  the  data. 

FSIG  - 

This  routine  calculates  the  variance  for  raw 
pseudorange  and  range  differenced  data  for  the 
satellite  being  tracked. 

GDOP  - 

This  routine  if  turned  on,  writes  a  GDOP  file. 

LINV3F  - 

This  routine  is  in  the  IMSL  and  performs  a 
matrix  inversion 

KALFIL  - 

This  routine  uses  corrected  data  to  do  the 
Kalman  filter. 

INKAL  - 

This  routine  sets  up  initial  covariance  matrix 
in  UDU'  format. 

UDFACT - 

This  routine  puts  Apriori  covariance  matrix 
into  UDU'  factorized  matrix. 

RKAL  - 

This  routine  calculates  partial  derivatives  of 
the  data  equation  for  pseudorange  data. 

UPDATE  - 

This  routine  will  update  the  Kalman  solution 
for  one  data  point. 

DPKAL  - 

This  routine  calculates  partial  derivatives  of 
the  data  equation  for  delta  (Doppler)  range 
differences. 

CONSTR  - 

This  routine  will  constrain  the  solution  to 
the  surface  of'  a  sphere  (height)  and  constrain 
velocity. 

XYZLLII  - 

Converts  X,  Y,  Z  in  an  earth-fixed  geocentric 
coordinate  system  to  geodetic  coordinate  system. 

PRNSOL  - 

This  routine  prints  solutions. 

COVEVL - 

This  routine  finds  the  covariance  matrix  from 
the  UDU'  factorized  matrix  in  U. 

BLOCK  DATA 

-  This  routine  was  added  at  NPS  to  handle 
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data  statements  and  common  statements  that 
had  ‘to  be  moved  in  the  code  or  written,  to  ^ . 
be  compatible  with  the  IBM-3033  FORTRAN  77 
compilers. 


MAIN  PROGRAM  KALMN2 
CARDIO 

LLHXYZ 

SIGXYZ 

MATRIX 

LINV2F 

VMULFF 

VMULFM 

INDAT 

TIMEUP 

PROMPT 

GETDAT 

DEDPDT 

DEFBR 

CRD  AT 

BSADR 

CRDAT2 

NEWRAN 

10NRD 

ION 

IONRD 

EPIIDAT 

CLOCK2 

CLOCK 

SATPOS 

ESTPOS 

ZANG 

RES 

SMOT 

TROP 

REL 

FSIG 

GDOP 

LINV3F 

KALFIL 

INKAL 

UDFACT 

LLHXYZ 

RGKAL 

UPDATE 
DPKA  L 

UPDATE 
CONST  R 

UPDATE 

XYZLLH 

PRNSOL 

COVEVL 

BLOCK  DATA 
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APPENDIX  G 

INPUT  GUIDE  FOR  PROGRAM  KALMN2 


(  by  S.  L.  Meyerhoff,  NSWC) 

Var  Type 

Description 

CARD  ONE 

col  1  to  10 

Format  ( 2F10.0.6I5 ) 

(This  card  can  be  read  format  free) 

TSS 

real  The  GPS  time-tag  for  data  at  start 

of  solution  (seconds) 

col  11  to  20 

TES 

real  The  GPS  time-tag  of  last  data  point 

to  use  in  solution  (seconds) 

col  21  to  25  (right  justified) 

NPN  integer  The  number  of  time  lines  between 

solutions 


col  26  to  30  (right  justified) 

NPB  integer  The  number  of  solutions  between 

solutions  being  printed 

col  31  to  35  (right  justified) 

ITRY  integer  The  number  of  solutions  printed 

covariance  matrix  prints 


col  36  to  40  (right  justified) 

NSAT  integer  The  number  of  Packers  to  process 

in  the  solution  (range  I  to  4) 


col  41  to  45  (right  justified) 

IG1  integer  =  1  if  data  source  is  standard 

9-track  format 
=  2  if  data  source  is  other 


col  46  to  50  (right  justified) 
ISP  integer 


col  51  to  55  (right  justified) 
ICLOCK  integer 


col  56  to  60  (right  justified) 
1AGE  integer 


=  1  for  satellite  position  from 
NAVIGATION  DATA  MESSAGE 
=  2  for  satellite  position  from 
Celestial  Earth  Fixed 
Trajectory  files 


=  0  clock  corrections  are  from 
NAVIGATION  DATA  MESSAGE 
=  1  clock  corrections  arc  post 
fit  data  from  cards 


=  1  make  GDOP  file 
otherwise  no  file  is  made 


col  61  to  65  (nght  justified) 

IAGE  integer  =  1  solve  for  clock  aging 

otherwise  no  file  is  made 


CARD  TWO  ( Receiver  Position)  Format  (/ W,7G  10.0) 

(This  card  can  be  read  format  free) 


col  10 
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I  FLAG 


integer 


=  0  for  Cartesian  coordinates 
=  1  for  Geodetic  coordinates 


col  11  to  20 

XL  AT 

real 

latitude  or  X  coordinates  for 
position  (degrees  or  kilometers) 

col  21  to  30 
YLOG 

real 

longitude  or  Y  coordinate  for 
position  (degrees  or  kilometers) 

col  31  to  40 
ZHT 

real 

height  or  Z  coordinate  for 
position  (kilometers) 

col  41  to  50 

TB 

real 

time  bias  for  receiver  clock 
(kilometers) 

col  51  to  60 

TD 

real 

time  drift  for  receiver  clock 
(kilometers/second) 

col  61  to  70 

TO 

real 

epoch  for  clock  drift  (given  and 
solution)  if  negative,  then  first 
used  time  line  is  the  epoch  for 
the  solution  clock  drift 

col  71  to  80 
AGE 

real 

aging  for  receiver  clock 
(kilometers/second"second) 

CARD  THREE  {Reference  Geo  id)  Format  {2G  10.0) 

(This  card  can  be  read  format  free) 

col  1  to  10 

AC  real  semi-major  axis  of  the  reference 

geoid  (kilometers) 

col  11  to  20 
ROBL 

real 

obliquity  of  the  reference  geoid 

CARD  FOUR  ( Sigmas  for  Position )  Formal  ( 7G10.0 ) 

(1  his  card  can  be  read  format  free) 

col  1  to  10 
SIGLT 

real 

sigma  of  X  or  latitude  coordinate 
(kilometer) 

col  11  to  20 
SIGLG 

real 

siema  of  Y  or  longitude  coordinate 
(kilometer) 

col  21  to  30 
SIGHT 

real 

sigma  of  Z  or  height  coordinate 
(kilometers) 

col  31  to  40 
SIG(i) 

real 

sigma  on  refraction  scaling  constant 

col  41  to  50 
SIG(2) 

real 

siema  on  receiver  time  bias 
(kilometer) 

col  51  to  60 
SIG(3) 

real 

sigma  on  receiver  clock  drift 
(kilometer/sccond) 

col  61  to  70 

SIG(4)  reai  sigma  of  velocity  in  X  direction 

(kiiomcters/second) 

CARD  FIVE  ( Sigmas  for  Velocity)  format  ( 7GI0.0 ) 

(This  card  can  be  read  format  free) 


col  1  to  10 
S1G(5) 

real 

sigma  of  velocity  in  Y  direction 
(kilometers/  second) 

col  11  to  20 
SIG(6) 

real 

sigma  of  velocity  in  Z  direction 
(kilometers/second) 

col  21  to  30 
SIG(7) 

real 

sigma  on  receiver  time  bias 
(kilometers) 

col  31  to  40 
SIG(8) 

real 

sisma  on  receiver  clock  drift 
(kilometers/second) 

col  41  to  50 
SIG(9) 

real 

sigma  not  used 

CARD  SIX  ( Types  of  solutions  to  do)  Format  (1615) 

(This  card  cannot  be  read  format  free) 

col  5 


NSOL(l) 

integer 

=  1  for  range  position  solution 
=  0  no 

col  10 

NSOL(2) 

integer 

=  1  Height  constraint  on  position 
=  0  constraint  only  if  less  than  4 
satellites 

col  15 

NSOL(3) 

integer 

=  1  write  position  solution  on  file 
=  0  no  solution  saved 

col  20 

NSOL(4) 

integer 

=  1  to.  solve  for  refraction 
correction  in  position 
solution 
=  0  no 

col  25 

NSOL(5) 

integer 

=  1  to  solve  for  refraction 
correction  in  velocity 
solution 
=  0  no 

col  30 

NSOL(6) 

integer 

=  1  to  solve  for  clock  drift 
from  position  solution 
=  0  no 

col  35 

NSOL(7) 

integer 

=  1  do  a  velocity  solution 
(this  has  not  been  used 
before) 

=  0  no 

col  40 

NSOL(8) 

integer 

=  1  do  a  time  update  of 
covariance  matrix 
=  0  no 
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col  45 

NSOL(9) 

integer 

=  1  not  used 
=  0  no 

col  50 

NSOL(IO) 

integer 

=  1  for  a  constraint  on 
velocity  solution 

=  0  no 

CARD  SEVEN  ( Variance  Bias  Factor)  Format  ( 2G10.0 ) 

(This  card  can  be  read  format  free) 

col  1  to  10 

K2 

real 

bias  factor  for  variance  on  range 
data 

col  11  to  20 

K4 

real 

bias  factor  for  variance  on 

Doppler  data 

CARD  EIGHT  ( Print  Options )  Format  (515) 

(This  card  cannot  be  read  format  free) 

col  5 

IPRT(l) 

integer 

=  1  print  position  solution 

(x,y;z)  f  . 

=  0  do  not  print 

col  10 

IPRT(2) 

integer 

=  1  print  solution  sigmas  and 
clock  terms 
=  0  do  not  print 

col  15 

IPRT(3) 

integer 

=  1  print  data  and  time  tags 
=  0  do  not  print 

col  20 

IPRT(4) 

integer 

=  I  residuals,  range,  trop, 
data  sigma 
=  0  do  not  print 

col  25 

IPRT(5) 

integer 

=  1  print  satellite  position 
=  0  do  not  print 

col  30 

IPRT(6) 

integer 

=  1  residuals  and  all  data 
corrections 
=  0  do  not  print 

col  35 

IPRT(7) 

integer 

=  1  print  position  in  lat,  long, 
and  height 

col  40 

IPRT(8) 

integer 

=  1  print  real  time  solutions  from 
9-track 

CARD  NINE  (weather  Datax  Elev.  Cutoff)  Format  (4G10.0) 

('1  his  card  can  be  read  format  free) 

col  1  to  10 

PB 

real 

pressure  in  millibars 

col  11  to  20 

TC 

real 

temperature  in  degrees  Centigrade 

col  21  to  30 


HD 

real 

relative  humidity  in  percent 

col  31  to  40 

EL 

real 

elevation  angle  cutoff  (degrees) 

CARD  TEN  ( Station  Number)  Format  (Free) 

(This  card  can  be  read  format  free) 

col  1  to  10 

Q(l) 

real 

time  propagation  noise  for  X 
velocity  (kilometers) 

col  11  to  20 

0(2) 

real 

time  propagation  noise  for  Y 
velocity  (kilometers) 

col  21  to  30 

Q(3) 

real 

time  propagation  noise  for  Z 
velocity  (kilometers) 

col  31  to  40 

Q(4) 

real 

time  propagation  noise  for  time 
drift 

col  41  to  50 
SIGPS 

real 

sigma  on  position  constraint 

col  51  to  60 
S1GVL 

real 

sigma  on  velocity  constraint 

CARD  ELEVEN  ( Diagonal  elements  of  state  transition  matrix) 

(This  card  can  be  read  format  free) 

col  1  to  10 

PH  1(1,1)  real  state  transition  matrix  for 

parameter  1 

col  11  to  20 
PHI(2,2) 

real 

state  transition  matrix  for 
parameter  2 

col  21  to  30 

PH  1(3,3) 

real 

state  transition  matrix  for 
parameter  3 

col  31  to  40 

PH  I  (4,4) 

real 

state  transition  matrix  for 
parameter  4 

col  41  to  50 
PH1(5,5) 

real 

state  transition  matrix  for 
parameter  5 

col  51  to  60 

PM  1(6,6) 

real 

state  transition  matrix  for 
parameter  6 

col  61  to  70 

PH  1(7,7) 

real 

state  transition  matrix  for 
parameter  7 

CARD  TWELVE  ( Diagonal  elements  of  state  transition  matrix) 

(This  card  can  be  read  format  free) 

col  1  to  10 

PH  1(8,8) 

real 

state  transition  matrix  for 
parameter  8 

col  11  to  20 
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PI-1 1(9,9) 

CARD  THIRTEEN 

col  1  to  10 
I-IT 


real  state  transition  matrix  for 

parameter  9 

(Height  for  position  constraint ) 

('Inis  card  can  be  read  format  free) 


real  Height  above  reference 

ellipsoid  to  constrain 
position  solution  to 
(kilometer) 
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APPENDIX  H 

ABBREVIATIONS  AND  ACRONYMS 


GPS 

— 

Global  Positioning  System 

NPS 

= 

Naval  Postgraduate  School  Monterey  CA 

DMA 

= 

Defense  Mapping  Agency 

NO  A  A 

National  Oceanic  and  Atmospheric  Administration 

NSWC 

= 

Naval  Surface  Weapons  Center 

SV 

= 

Space  Vehicle 

SEP 

= 

Spherical  Error  Probability 

P-code 

= 

Precision  code 

C/A-code 

= 

Coarse  Acquisition  code 

WGS 

= 

World  Geodetic  System 

PDOP 

= 

Position  Dilution  of  Precision 

HDOP 

= 

Horizontal  Dilution  of  Precision 

VDOP 

= 

Vertical  Dilution  of  Precision 

TDOP 

= 

Time  Dilution  of  Precision 

SEF 

=r 

Standardized  Exchange  Format 

CIO 

= 

Conventional  International  Origin 

ARL 

= 

Applied  Research  Laboratories  at  the  University  of  Texas 
at  Austin 

IBM 

= 

International  Business  Machines 

CDC 

= 

Computer  Data  Corporation 

m 

as 

meters 

km 

= 

kilometers 

s 

= 

second 

Hz 

= 

Hertz 

MHz 

= 

Megahertz 
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